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Preface 


This  report  is  the  result  of  my  investigation 
into  numerical  solutions  of  differential  equations, 
expressed  in  both  differential  and  integral  form. 

The  method  of  finite  differences  was  utilized  for 
solutions  of  the  equations  in  differential  form.  A 
Fredholm  integral  representation  was  developed  and 
numerically  evaluated  to  compare  with  the  finite 
difference  method.  The  purpose  of  this  thesis  was 
to  examine  the  solutions  obtained  to  the  differen¬ 
tial  equation  using  both  above  methods  and  report  on 
the  advantages  and  disadvantages  of  each.  A  secondary 
purpose  of  this  thesis  was  to  strengthen  my  own  mathe¬ 
matical  background  on  numerical  methods  and  their 
ability  to  solve  differential  equations.  I  have 
attempted  to  include  sufficient  detail  to  provide  the 
reader  with  a  step  by  step  account  of  the  development. 

I  want  to  thank  Dr.  Bernard  Kaplan,  my  advisor, 
for  his  guidance  and  assistance  throughout  this  study. 

I  would  also  like  to  thank  Dr.  John  Jones  for  his 
many  helpful  discussions  pertaining  to  spline  theory 
and  also  to  Dr.  Wilhelm  Ericksen  for  his  patience  in 
helping  correct  many  of  my  computer  programs.  I  also 
wish  to  thank  Dr.  W.  Kessler  of  the  Air  Force  Materials 
Laboratory  for  sponsoring  this  study.  Finally,  I 
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Abstract 


A  study  of  several  numerical  methods  for  the 
solution  of  boundary  value  problems,  in  both  one  and 
two-dimensions,  was  conducted  using  the  CDC  66 00  com¬ 
puter.  The  method  of  finite  differences  was  employed 
for  solution  of  the  equations  in  differential  form. 

These  numerical  solutions  were  compared  to  those  ob¬ 
tained  by  transforming  the  original  differential  equa¬ 
tion  into  integral  form  and  approximating  their  solu¬ 
tion  using  numerical  integration  via  the  trapezoid 
rule.  All  numerical  experiments  were  conducted  using 
Dirichlet  boundary  conditions.-^ 

w^In  the  one-dimensional  cases  studied  it  was  found 
that  both  methods  are  equivalent,  i.e.,  yield  identi¬ 
cal  solutions  wheft  the  integral  representation  had  a 
linear  weighted  Green's  function  kernel.  In  addition, 
the  integral  approach  was  found  to  be  as ) accurate  in 
all  one-dimensional  cases  as  the  mettxfd  of  finite 
differences.  The  finite  difference  method  proved 
quicker  than  the  numerical  integration  techniques  in 
all  but  one  test  case  where  the  Green's  integral  rep¬ 
resentation  was  examined. 

For  the  two-dimensional  investigation  the  steady- 
state  heat  conduction  equation  was  analyzed.  Again, 
the  method  of  finite  differences  in  two-dimensions _ ^ 


x 


was  compared  to  the  integral  approach,  using  cubic 
splines.  The  method  of  finite  differences  was  found 
to  be  superior  in  calculating  the  internal  temper¬ 
ature,  at  all  nodal  points,  as  compared  to  the  inte¬ 
gral-spline  solution. 


I .  Introduction 


Background 

A  majority  of  problems  encountered  in  techni¬ 
cal  research  by  engineers  and  physicists  can  be  ex¬ 
pressed  in  mathematical  form  as  a  differential  or  par¬ 
tial  differential  equation.  Solutions  of  these  equa¬ 
tions  are  dependent  upon  initial  conditions  and/or 
boundary  values.  (If  the  values  of  the  function  are 
specified  on  the  boundary,  the  equation  is  said  to 
contain  Dirichlet  boundary  conditions.  If  the  normal 
derivatives  (gradients)  of  the  function  are  specified 
on  the  boundary  then  Neumann  boundary  conditions  are 
said  to  exist.  The  boundary  conditions  are  referred 
to  as  mixed  if  the  initial  conditions  describing  the 
differential  or  partial  differential  equation  contain 
both  Dirichlet  and  Neumann  conditions.)  Because  the 
equations  possess  unique  solutions  it  does  not  nece¬ 
ssarily  follow  that  these  solutions  are  easy  to  ob¬ 
tain.  In  many  cases  the  exact  closed-form  solutions 
are  not  attainable;  and  as  a  result,  approximation 
techniques  must  be  used  to  generate  analytical  values. 
For  this  reason,  most  of  all  important  problems  re¬ 
quire  application  of  some  numerical  method. 

Numerical  Techniques 


Though  there  are  a  variety  of  techniques  available 


|  to  handle  specific  boundary  value  problems,  this 

thesis  will  be  concerned  with  three  specific  methods: 

(1)  The  method  of  finite  differences 

(2)  Fredholm  integral  equations  and  their  numeri¬ 
cal  approximations 

(3)  utilization  of  cubic  splines 

One-Dimensional  Cases 

The  method  of  finite  differences  is  based  upon  a 
scheme  of  numerical  differentiation.  The  original 
differential  equation  is  replaced  by  a  finite  number 
of  algebraic  expressions  defined  over  the  interval; 
this  set  of  equations  is  easily  solved  using  the  com¬ 
puter.  In  the  one-dimensional  cases  methods  (1)  and 
(2)  are  compared  for  both  accuracy  and  speed  of  com¬ 
putations. 

The  second  method  transforms  the  original  differ¬ 
ential  expression  into  an  equivalent  integral  equation, 
usually  a  Fredholm  integral  equation  of  the  second  or 
third  kind  (Ref  1:  381-382).  The  advantages  of  this 
technique  is  that  in  many  cases  a  weighted  symmetric 
Green's  function  kernel  results  and  can  be  evaluated 
using  standard  numerical  integration  techniques,  such 
as  the  trapezoid  rule. 

( 
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Two-Dimensional  Cases 


In  the  two-dimensional  analysis  the  steady- state 
heat  conduction  problem,  with  an  inhomogeneous  boun¬ 
dary  condition  of  along  one  edge,  was  investigated. 
The  analytical  solution  was  calculated  and  compared 
at  16  interior  nodal  points  with  the  method  of  finite 
differences . 

The  third  method  investigated  in  the  two-dimen¬ 
sional  case  was  again  one  of  converting  the  original 
partial  differential  equation  into  integral  form  by  a 
method  proposed  by  Hajdin  and  Krajcinovic  (Ref  2t  523- 
539)*  The  unknowns  appearing  within  the  integrals 
are  approximated  by  cubic  splines  and  numerically 
integrated.  Again  methods  (1)  and  (3)  were  compared 
for  both  accuracy  and  ease  of  computation. 

Purpose 

There  are  several  important  for  investigating 
numerical  integration  techniques  two  methods  bases  upon 
numerical  differentiation.  Computer  utilization  costs 
are  normally  are  direct  function  of  utilization  time. 
Thus  it  becomes  of  prime  importance  to  use  the  most 
cost  effective  method  in  solving  a  particular  problem. 
Hajdin  and  Krajcinovic  contend  that  numerical  inte¬ 
gration  is  many  times  more  accurate  than  numerical 
differentiation  and  that  results  through  numerical 
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integration  should  be  superior  to  those  based  upon 
the  method  of  finite  differences.  For  the  same  accu¬ 
racy  one  should  be  able  to  significantly  reduce  the 
number  of  points  (i.e.  number  of  algebraic  equations) 
with  numerical  integration  (Ref  3*  509-510)*  The 
purpose  of  this  study  is  therefore  to  determine  whe¬ 
ther  or  not  numerical  integration  is  advantageous  or 
comparable  to  the  method  of  finite  differences,  a 
method  of  numerical  differentiation. 


Plan  of  Development 

Due  to  time  restrictions  placed  on  this  study, 
four  types  of  problems  were  considered.  In  the  one¬ 
dimensional  analysis  the  problems  and  boundary  con¬ 
ditions  were  the  following* 

(1)  y"  +  y  =  o,  y(o)»o,  y(»)»  \ 

(2)  if  +  =  1 »  ^/C°)=  o, y •)  =  o 

(3)  =  **J  yo*o,ys>*o 

Equation  (1),  equivalent  to  the  one-dimensional  Helm¬ 
holtz  equation,  is  introduced  to  illustrate  the  pro¬ 
cedures  used  in  transforming  differential  equations 
into  their  equivalent  integral  representation.  In 
addition,  equation  (1)  is  also  used  to  demonstrate 
the  numerical  methods  employed  to  evaluate  the  Fredholm 
integral  representations  of  equations  (2)  and  (3). 


All  equations  had  closed-form  analytical  solu¬ 
tions  that  were  used  for  comparison  to  the  numerically 
generated  values.  Numerical  computation  was  carried 
out  on  the  CDC  6600  computer  using  identical  numeri¬ 
cal  techniques  in  all  cases.  Also,  the  number  of 
iterations  per  interval  was  kept  the  same  so  that  an 
accurate  comparison  of  both  methods  could  be  made. 
Solutions  were  first  obtained  for  1-10  iterations 
over  the  boundary,  and  increased  to  50  iterations/in¬ 
terval  for  the  final  analysis.  All  computation  was 
carried  to  six  significant  decimal  places. 

In  two-dimensional  analysis,  Poisson's  equation 
was  investigated.  The  problem  chosen  was  the  steady- 
state  heat  conduction  over  a  square  plate  with  an 
inhomogeneous  Dirichlet  boundary  condition  at  one  edge. 
The  known  analytical  solutions  were  developed  using 
the  method  of  separation  of  variables.  Over  10,000 
independent  series  summations  were  required  at  each 
interior  node  point  to  yield  six  decimal  place  conver¬ 
gence.  The  method  of  finite  differences  was  compared 
with  the  analytical  results  at  16  and  81  interior 
points.  For  integral  conversion,  the  method  of  cubic 
splines  was  used  and  numerically  integrated.  The  cal¬ 
culated  solution  was  compared  with  the  known  analytic 
solution. 
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II.  Theory 


Integral  Equations 

An  integral  equation  is  one  in  which  the  function 
to  be  determined  appears  under  an  integral  sign. 

Linear  integral  equations,  that  is,  equations  in 
which  the  unknown  function  -f  appears  to  no  higher 
power  than  one,  are  conventionally  divided  into  two 
classifications.  An  equation  of  the  form 

<*(x)f(x)  =  fU)  +  XJ  K(*.$)f(C)  (2) 

a. 

where  <K,  Jr,  and  K  are  given  functions  and  and 

b  are  constants  is  known  as  a  Fredholm  equation. 

The  function  K(x,s)  is  known  as  the  kernel  of  the 
integral  equation  and  is  frequently  a  weighted  Green’s 
function.  If  the  upper  limit  of  the  integral  is  not 
a  constant  the  equation  takes  the  form 

<*(x)-f<x)  =  Foe)  +  ^  \  K<x»^>“fwds  (3) 

J<K 

and  is  known  as  a  Volterra  equation. 

When  ^  ^  0  *  equation  (2)  involves  the  unknown 
function  both  inside  and  outside  the  integral.  If 
0(*O,  the  unknown  function  appears  only  under  the  in¬ 
tegral  and  the  equation  is  known  as  an  integral  equa¬ 
tion  of  the  first  kind.  If  PC  =  l  the  equation  is 
said  to  be  of  the  second  kind.  In  the  more  general 
case  whenaf(X}  is  not  a  constant  the  equation  is 
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called  an  integral  equation  of  the  third  kind  (Ref  It 

382). 


Multiple  Integrations 

Certain  integral  equations  can  be  deduced  from  ' 
or  reduced  to  differential  equations.  In  order  to 
accomplish  the  reduction  it  is  frequently  necessary 
to  make  use  of  the  formula* 


x  x  x  A  J* 

J TfCxjdx^X  r  Qx J'Tf'njdn  = 

CL  ft  ft  CL 

Equation  (4)  is  obtained  by  integrating  the  left-hand 


side  by  parts,  that  is 

f*x 


f  f=*  X  $  A  * 

£§jV<i)dn7  -  Jffjj: M 

*  a  a  4 

More  generally,  by  a  repeated  application  of  this 

procedure,  the  results  of  integrating -ft*)  n  times 
over  the  same  limits  is 


(5) 


n  iv»es 

Expression  (6)  will  be  useful  in  manipulating  multiple 
integrations  in  the  work  which  follows  (Ref  4*  722). 


Differential  Conversion  to  Integral  Form 

To  illustrate  the  mechanics  for  converting  a 
differential  equation  into  integral  form,  consider 
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the  boundary  value  problem 

p£  =  -foo 

dX4 

with  boundary  conditions 


=  a 


&)*  b 


(7) 


(8) 


First  integrate  both  sides  of  (  7  )  with  respect  to 
X  over  the  interval  (0,x) 


or 


b. 

<Jx  i o 


o 

X 


=  -  ffcx)d 

o  J 


(9) 


(10) 


Therefore  (  7  )  is  equivalent  to  the  expression 


cjjpo  -  c  -  fa*)** 


(11) 


where C  represents  the  unknown  value  of 
A  second  integration  over  ( 0,X  )  leads  to 

=  Cfdx  -  J  j'fcxjdxdx 


or 


X.X 


o  O  o 


J; 


X  dX 


(12) 


(13) 


This  becomes 

|X 


•y 

«  i  i 

Using  equation  (  6 )  for  the  right-hand  integral  and 


^<*)|  =  Cx  ~  |  f  f(x)*X*X 


(14) 
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simplifying,  equation  (14)  becomes 

j(x)=  y(o)  +  Cx  (15) 

or 

X 

wj(x)  ft  4  Cx  -  )  <4  d6) 

The  constant  C  can  be  evaluated  by  applying  the  sec¬ 
ond  boundary  condition,  namely 

g 

y(t)s  \>  = «. ■*  cJt  -  ut-eM«>a*  a?) 

Solving  for  C  gives 

,  ,  ri 

c  =  ^ 4  jr  j  d8) 

Substituting  C  above,  equation  (16)  becomes 

.J  ^ 

yCx)*  a  <  19) 

This  last  expression  can  be  rewritten  by  expanding 

X  f 

the  first  integral,  i.e.,J=j  +J  .  Equation  (19)  then 

,  o  o  x 

becomes 


*J(x)=  a  +  j.  j(J  -fK(f)  is 

0 

r*  J 


r  • 


Collecting  terms  under  the  same  integral  limits,  (20) 


becomes 


r  x 


+  «• 


(21) 


With  the  abbreviation 


K(*i«) 


Jfd-f) 

j(l-X) 


equation  (21)  becomes 


*>? 


(22) 


yw  s  J  K<x.$MVO  <1$  +  4-  a  (23) 

Thus,  the  integral  equation  corresponding  to  the 
boundary  value  problem  (7)  is  a  Fredholm  equation  of 
the  second  kind. 


Kernel  Properties 

Note  that  the  kernel,  K(x.O  given  by  equation 
(22),  has  different  analytical  expressions  in  the 
two  regions,  X<  $  and  X>£  ,  but  that  the  expressions 
are  equivalent  when  X-f  .  Observe  also  that  in  each 
region  K  is  a  linear  function  of  X  and  that  If  van¬ 
ishes  at  the  end  points.  K(*,n  can  be  thought  of  as  a 
function  of  X  for  a  fixed  value  of  £  . 
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Fig.  1  The  Integral  Kernel,  K(x,s) 


Finally,  K(x,$)  is  unchanged  if  X  and  %  are  inter¬ 
changed;  that  is,  =  Ktex)  •  Kernels  having 

this  last  property  are  said  to  be  symmetric  and  obey 
the  Reciprocity  principle.  Many  of  the  kernels  gen¬ 
erated  in  converting  boundary  value  problems  into  in¬ 
tegral  equations  have  this  property  and  are  commonly 
identified  as  Green's  functions. 

The  Green's  Function 

Green's  functions  are  extremely  important  mathe¬ 
matical  quantities  in  physics.  Mathematically,  they 
allow  the  solution  of  differential  equations  to  be  ex¬ 
pressed  in  integral  form.  The  solution  of  some  inhomo¬ 
geneous  differential  equation 

L  ‘j(x)  t  f(x)  =  o  (24> 

with  homogeneous  boundary  conditions 


11 


(2J) 


y(o)=0,  MCJ1)=0 

and  where  L_  is  the  self-adjoint  differential  op¬ 
erator, 

L=  + 

can  he  expressed  in  the  integral  form 

A 

y(x)=  j  G(x,?)f(5)  (27) 

o 

where  G(*,5)  ■is  called  the  Green's  function  (Ref  6:  599- 
600) .  The  advantages  of  the  Green's  function  are: 

(1)  The  homogeneous  boundary  conditions  are  in¬ 
corporated  in  the  integral  representation  (27). 

(2)  Equation  (2?)  enables  the  differential  rep¬ 
resentation,  (24),  to  be  solved  via  the  inte¬ 
grating  process. 

(3)  The  Green's  function  method  also  allows  solu¬ 
tions  of  (24)  with  inhomogeneous  boundary 
conditions  (Ref  5*  12-14). 

Determining  the  Green's  Function 

The  Green's  function  satisfies  the  differential 
operator  of  a  homogeneous  boundary  value  problem, 
everywhere  except  at  one  point.  In  addition,  it 
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jrv  i  • 


vanishes  at  the  end  points  and  has  a  discontinuous 
first  derivative  at  X  -  %  •  For  a  boundary  value  prob¬ 
lem  of  the  form 


L  ’J(x)  r  -  -f(x)  <28) 

the  Green's  function,  Gfx.O  ,  is  defined  as  the  solu¬ 
tion  of  the  equation 


L  G(x,5)  =  -$<*-«) 


(29) 


where  £(*-?)  is  known  as  the  Dirac  delta  function 
having  the  properties 


r 


o 

CO 


x  -  f 


(30) 


The  Green's  function  satisfies  the  homogeneous  diff¬ 
erential  operator,  L.  ,  at  all  other  points  other  than 
.  At  X=5  a  singularity  exists  and  is  governed  by 
the  properties  of  the  delta  function  (Ref  15s  7). 

Over  the  interval  (0,fi)  it  is  possible  to  obtain  a 
convenient  form  of  a  solution  for  G(x,f)by  assuming 
that  the  Green's  function  can  be  represented  by  G,(x) 
whenX<$»  and,  by  Gj(X)  when  X>§  •  In  the  one-dimen¬ 
sional  case  the  Green's  function  has  the  following 
properties! 
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Cg ,<*) 

(1)  GM  =  ]  •  and  the  functions  G,(x)  and 

£G,uq  x>i 

G *(*)  each  satisfy  equation  (29)  in  their  in¬ 
tervals  of  definition,  i.e., 


LG, 00  =  O-nd  L  G2(x)=  0  (-for  X*f)  (3D 

(2)  The  function  G(x»0  satisfies  homogeneous  con¬ 
ditions  at  the  interval  end  points  o  and  X=Jl  . 

G,(°)=o,  6SU)  =  0  <32) 

(3)  G(x,0  is  continuous  at  X  =  "^  .  This  requires 
that 


G,($)  =  G^s) 

(4)  The  derivative  of  Gfof)has  a  discontinuity 
of  magnitude  -  tyco  at  Xs §  »  that  is, 


^<X)|  -  =  -  y, 

x-*5  4x  x-*S  * 


f<f) 


5 

Using  these  four  conditions,  it  is  possible  to  deter¬ 
mine  G(X,$)  and  represent  the  solution  of  the  differ¬ 
ential  equation  in  integral  form. 


(33) 


(34) 


Example 

As  an  example  of  determining  the  Green's  function 
for  a  differential  operator,  consider  the  boundary 
value  problem 
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f 


(35) 


with  inhomogeneous  boundary  conditions 


( o)  -  =  b 


(36) 


The  Green's  function  is  determined  as  the  solution  of 


L  &<*<«)=  <37> 

where)  si-  in  this  case.  Equation  (37  )  becomes 

££<*,§)  =  -  C(x-o  (38) 

<Jx*  d  v 

Assuming  that  Gfof)  can  be  expressed  by  Qjx<x)  and  GtU) 
over  the  interval  (0,1)  ,  the  Green’s  function  will  be 
of  the  form 


G(x,?) 


Gjx)  x<$ 

G^x)  x>% 


For  £  .  6, GO  and  G2(x)  satisfy  the  homogeneous 
equation  (31  )  *. 


(39) 


L  G,(x)  =  L  6t(x)  =  O  (*0) 

The  simplest  form  of  solution,  satisfying  equation 
( 40  ) ,  is  to  assume  that  takes  the  form 
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~  l  A  x  +  B 

G^?)=  < 

Cx  +  D 


(41) 


X<c 
*>S 

From  the  previous  properties  (2),  (3)  and  (4)  for  the 
Green's  function,  all  of  the  constants  in  ( 41  )  can 
he  evaluated.  From  property  (2) ,  G(x,f)  must  vanish  at 
the  ends,  therefore 

G*(0)  r  O  =  &  (42) 

and 


G2U)  =  o  =  a  d=  -c*  (43) 

Substituting  these  values  for  b  and  [)  back  into  (41  ) 
the  Green's  function  becomes 


g<«i.  r*  *<( 

(44) 

[C(x-f)  x>§ 

Using  the  fact  that  G(x,f)is  continuous  at  X=  §  » 
perty  (3),  namely 

pro- 

G,U)  =  Gj?) 

(45) 

the  constant  Q,  can  be  expressed  in  terms  of  /]  s 

As  =c(f-s) 

(46) 

•  •  c . 

(47) 
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r 


Substituting  expression  (4?  )  back  into  (44  ),  the 


Green’s  function  takes  the  form 


G<x,f)= 


*L  U-o) 


(48) 


Using  property  (4),  the  fi:  s. 

be  discontinuous  at  X  =  f  1 

4S(X)|  “  d  G,(*)  1 

<ix  X-*f  <Ax  x-+f 

Doing  the  differentiation  yields 

M.  -  4  -/ 

or 

j.&zi) 

x 

and  equation  (48  )  becomes 


rivative  of  must 

fically, 


Equation  (52  )  represents  the  final  expression  for 
G(x,$)  for  the  operator  |_s  Notice  that  this  is  ex¬ 

actly  the  form  of  the  kernel  for  the  integral  expansion 
of  this  boundary  value  problem,  given  by  equation  (22  ). 


The  Green's  Integral  Equation 


Many  different  boundary  value  problems,  upon 
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integral  conversion,  have  symmetric  kernels  equivalent 
to  weighted  Green's  functions.  Therefore,  two  pro¬ 
cedures  to  generate  the  Green's  function  can  be  used: 

(1)  Through  integral  conversion  or 

(2)  Finding  G(x,f)  from  the  differential  operator,  L. 
It  should  also  be  mentioned  that  not  all  integral 
conversions  yield  symmetric  Green's  functions.  In 
particular,  if  the  differential  equation  is  of  the  form 

4^  0^3  +8  u  =  O  (53) 

with  homogeneous  boundary  conditions 

y(JO=o  (5^) 

then  the  resulting  kernel  of  the  integral  equation  is 
both  nonsymmetric  and  discontinuous  at  unless 

O  (Appendix  A). 

Both  methods  (1)  and  (2)  should  be  equivalent, 
however  different  expressions  forG^§)can  occur.  As 
an  example  consider  the  one-dimensional  Helmholtz 
boundary  value  problem 

4^  +  AN  =  I  <55> 

dx*  d 

with  boundary  conditions 

y(o)=  Oj  y<0=  o  (56) 

Upon  integral  conversion,  page  68  ,  it  was  found  that 
;/<*>  was  of  the  form 


I 

T)=  i)  <U 


(57) 


where  the  kernel  is  defined  by  the  expression 

[ x( i-f)  x  <5 

(58) 

l  SO-*)  x  > 5 

Expression  (  58)  represents  the  equivalent  weighted 
Green's  function  for  the  integral  conversion  method. 
Next,  consider  the  problem  in  operator  form,  where  L. 
represents  the  one-dimensional  Helmholtz  operator 

d1  '*• 


I  :f  +X 

L  dxl 


(59) 


The  Green's  function  must  obey  the  differential  equa¬ 
tion 

£^x,0  +  X’Gfx.O  =  -  $(x-f)  (6o) 

»Xl 

From  the  method  outlined  using  the  differential  operator, 
l_  ,  the  Green's  function  was  found  to  be  (Appendix  C)» 

*nXx  x<S 

G(x,o=V  x  u  «*> 

I  Sl<\  X(\-x)  „  V  _ 

and  M(*)  is  given  by 

J  r‘ 

y(X)  -  -  (62) 

where  6^f)  is  defined  by  equation  (61  )  and  f<5)  a - i  from 
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(5  5).  Clearly,  vj(x)  is  represented  equivalently  by  two 
completely  independent  forms  of  the  Green's  functions 
that  is,  equations  (57)  and  (62),  each  with  distinct, 
different  kernels,  must  be  equivalent  expressions. 


Equivalency  of  the  Integral  Equation 
and  the  Green's  Integral 

It  will  now  be  shown  that  equations  (57)  and  (62) 
are  identical  expressions  and  that  they  both  satisfy 
the  one-dimensional  Helmholtz  differential  equation. 
Equation  ( 57)  can  be  written  in  the  form 


o  J  'x 

where  the  kernel,  defined  by  equation  ( 58)  has  been 
substituted  over  the  appropriate  limits.  According 
to  Leibnitz's  rule  for  differentiation  under  an  integral, 
(Ref  li  383),  the  first  derivative  of  (63)  becomes 


,  1 

y<*)=  -hOfya-Odj  +  [(t-Ofryo-iMf 


(64) 


The  second  derivative  of  JW  is  found  by  differentiating 
(64)  with  respect  to  X  ,  and  is 


-X0tyx)-l )  -  (i-xK^yx)-i)  =  \  -  fytx)  (65) 

Substituting  the  equation  for  y'ix) ,  equation  (  65) ,  back 
into  the  original  differential  equation,  expression  ( 55) » 
it  is  found  that 
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y‘{x)+)ty*)5s  [  ^yx>J  *  >>  s 


Hence,  the  integral  representation  given  by  (57) 
satisfies  differential  equation  (55)  and  therefore 
must  be  equivalent  due  to  the  uniqueness  of  solution 
guaranteed  through  the  given  boundary  conditions. 

Equation  (62)  can  likewise  be  written  in  the  form 

X  ^ 

V*).  -  ic  -  faaMr-tf  »«»>«  J{  (6; 

o  X 

where  again  has  been  substituted  over  its  defined 
limits.  Using  Leibnitz's  rule  the  first  derivative 
of  y(x)  becomes 

tied-  A  f fcoS  (6) 

J  J  A**M  J  As,n>* 

o  X 

The  second  derivative,  differentiating  expression  (68), 


becomes  the  quantity 

\z  f»in  X(i-x)sin\f 

*  A  J^tai — 


5twAxsroAfi-f) 

A*»«AJI 


.  Xcq$MJ?-x)siis>Ax  -t  \cosAxsvnA(i-x) 

Asm  A  A 

Substituting  the  integral  expression  for  y<>»)  and  vj(x) 
back  into  the  original  differential  equation,  it  is 
found  that 

4  0c)*A\<*)-  A  COS  Mt-  X)  Ax  +•  >  cos  A  x  s m  A(l-x) 

J  J  A  sm  XA 
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After  simplifying,  the  right-hand  side  of  equation  (70) 
becomes 

m'«h>v>=  Al"  z  i  (n) 

J  3  Xs>«X* 

Therefore,  the  Green’s  integral  representation  also 
satisfies  differential  equation  (55)  and  both  integral 
forms,  equations  (57)  and  (62)  must  be  equivalent  ex¬ 
pressions  . 

Integral  Representation 

One  last  procedure  will  be  introduced  to  demon¬ 
strate  how  equation  (62)  originates.  Consider  the  two 
equations 

(72) 


ja 

s  --f(x)  ,  y(o)*at  y(Z)=  b 


and 


iz G(x,s)  _ 


(73) 


Equation  (73)  represents  the  Green’s  function  solution 
for  the  differential  operator  L5  •  Multiply  equa¬ 
tion  (72  )  by  GOc,$)  and  equation  (73  )  by  y<x)  ,  subtract 
and  integrate  over  the  interval  (CJJ?)  to  obtain 


r»  i 


(7^) 
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Rewriting  the  left-hand  derivatives  in  equation  (7*0 
and  using  the  integral  properties  of  the  Dirac  delta 
function  (Ref  5*  6),  equation  (7*0  becomes 


'  vklW**)}**  = 

”  fG(x,0  fuxjix  +  yro 


(75) 


Integrating  the  left  side  by  parts,  using  the  fact  that 

JuJv=  =  uv-^vJadx  (76) 

The  left-hand  side  of  equation  (75)  becomes 

i <• J0A 

Evaluating  at  the  limits,  (77)  becomes 


(77) 


GttOjM)  -  6(0,1)  ^  -  yijc&ixt)  \  +  |  (78) 

Because  the  Green's  function  is  zero  at  end  points  the 
first  two  quantities  on  the  left  drop  out.  The  only 
expressions  left  to  evaluate  are  the  derivatives  of 

on  the  appropriate  intervals.  Differentiating 
equation  (52),  which  is  the  known  G(Kf)  for  the  operator 
fe.  .  gives 


Ir  t  X  <  § 

dG  (*,f)  -  J  5 

dx  J  -f 

t  r  *>$ 

Using  cJ6fy,f)  above,  equation  (78)  becomes 


(79) 
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(80) 


(81) 


Solving  for  ^(|)»  and  substituting  the  boundary  con¬ 
ditions  at  yo)  and  V*)  ,  equation  (80  )  becomes  the 
expression 

r* 

y(§)=  \  G<*f)'f(ioJx  +  -V- 

0 

Interchanging  the  labels  x  and  |  from  the  Reciprocity 
principle  (Ref  6s  328)  and  using  the  symmetry  of  G<x,f), 
equation  ( 81 )  is  transformed  into 

rJ 

J  G<*iS) ■?<%)«!§  +  k*  +  (82) 

This  equation  is  the  Green's  integral  solution  of 
differential  equation  (72),  with  inhomogeneous  Dirichlet 
boundary  conditions.  Note  the  similarity  of  this  rep¬ 
resentation  versus  the  integral  equation  expression 
given  by  equation  (23). 


Numerical  Integration 


b 


An  integral  of  the  form  stands  for  the 

area  represented  by  fix)  over  the  interval  («,*).  This 
area  can  be  approximated  by  several  methods;  the  tra¬ 
pezoid  rule  is  based  on  linear  interpolation,  the 
integral  is  expressed  as  a  sum  of  N  trapezoidal  areas; 
Simpson's  rule  for  numerical  integration  is  based  upon 
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quadratic  interpolation,  parabolas  or  polynominals  of 
degree  2  or  less  are  used  to  approximate  Tex)  between 
three  successive  points;  the  Newton-Cotes  method  is 
based  upon  polynominal  interpolation  of  degree  3  or 
more.  In  this  thesis  the  trapezoid  rule  is  used  to 
numerically  integrate  all  one-dimensional  cases. 


Trapezoid  Rule 

The  area  of  a  typical  trapezoid  with  base  length 
\*\  and  sides  fto;.,)  .  and  f<*  ;)  is 

/>rea  --  £[f(x^)+  (83) 

The  combined  area  of  all  trapezoids,  over  the  interval 

0.=  x,  to  t>?X^is 


T«*al 


area 


=  h 


+  (84) 
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Therefore,  for  the  trapezoid  rule  equation  (84)  becomes 

J'fooAx^  h£f<x;)+24Vxo+  •  •  •  +2f(xvj_,)  +  f (85) 

<L 

where  X,  ,  X2 . X^  are  equally  spaced  points  dividing 

the  interval  from  4*  x,  to  \>a  into  (N-l)  equal  parts, 
of  length  h  =  The  trapezoid  rule  can  be  made 

as  accurate  as  desired  by  choosing  h  very  small,  i.e., 

N  very  large.  In  addition,  there  are  several  advantages 
in  using  trapezoid  interpolation  versus  the  Simpson's 
rule. 


Trapezoid  vs.  Simpson's  Rule 

It  was  shown  earlier  that  all  differential  equa¬ 
tions  of  the  form 


ST  r  "  f  (x> 


with  homogeneous  end  conditions  have  linear  Green's 
kernels  when  expressed  in  integral  form  (Ref  li  461- 
462).  The  linear  kernel  of  equation  (86  )  is  shown  in 
Fig.  1  and  given  by  equation  (22).  Due  to  the  linearity 
of  K(M)  and  the  presence  of  a  'corner'  at  X  =  £  ,  Simpson’s 
rule  for  numerical  integration  will  be  less  accurate 
than  the  trapezoid  rule  (Ref  7*  217,  Problem  #4-25). 
Trapezoids  furnish  better  approximations  to  linear 
kernels  than  parabolas.  Parabolas  always  exceed  the 
actual  area  bounded  by  three  consecutive  points  across 


r 


\ 


the  discontinuity,  at  X=f  ,  for  a  linear  kernel.  In 


addition,  Simpson's  rule  can  only  be  used  if  N,  the 
number  of  points  in  the  interval  (a,k),  is  odd.  By 
using  Simpson's  rule  with  N  even  the  computed  integral 
will  exceed  the  true  value  because  of  the  multiple 
contributions  from  trapezoids  adjacent  to  the  dis¬ 
continuity.  As  an  example  to  illustrate  what  happens 
by  using  Simpson's  rule  when  N  is  even,  consider  the 
area  bounded  by  the  linear  kernel  given  by  equation 
(22)  over  the  interval  (0,1).  The  linear  kernel  is 

Uv  CX('~S)  X<S 

K(X,£)  =  <  (87) 

%  (  »  -  X)  *  >£ 

Divide  the  interval  (0,1)  into  three  regions  defined 
by  the  four  equally  spaced  points 

vb,  =  v- ' 

Fixing  and  applying  Simpson's  rule  to  the  area 

defined  by  one  extra  contribution  from  the 

region  defined  over  (Xj^Xj)  is  added  into  the  total  area 


Fig.  3  The  Linear  Kernel,  he*,*) 


The  first  contribution  comes  from  applying  quadratic 
interpolation  over  the  first  three  points  (XwXa/x3  ). 
Recall  that  in  quadratic  interpolation  a  polynominal 
of  degree  2  or  less  is  used  to  represent  a  function  "fix ) 
over  two  successive  intervals  (X„Xj,Xj).  Mathematically, 
if  the  value  of  fix)  is  known  at  three  locations,  a 
polynominal,  f<*?  ,  can  be  determined  of  the  form 

p(x)  =  A  X1  ■*  Bx  *-  C  (88) 

/ 

where  6,  and  C  are  all  constants  to  be  determined 
from  the  three  points 

(X'J(*'))}  (xiA<*t)]}  fa.ftxj)) 

Once  determined,  f6c)  is  then  substituted  for  'ffx)  over 
the  interval  (Xi,  X»,  Xj  )  and  integrated  by  expression 
(85).  Over  the  points  X„  X{,  and  given  above,  and 
using  the  value  of  on  the  interval  X<$,  f>(*)can  be 

found  and  is 

p(x)  a  2L  (89) 

Therefore,  Simpson's  method  applied  to  the  first  three 
points  in  Pig.  3  is  the  area  represented  under  the  line 

K6*f>  V3  . 

To  calculate  the  remaining  area  under  K(x,f)for  X>J  , 
three  additional  points  and  the  values  of  must  be 
known.  Now  however,  only  two  points  remain  for  the  in- 
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terval  x>f  ,  Xj  and  .  If  x4  is  used  again,  in  con¬ 
junction  with  Xj  ,  and  ,  a  parabola  can  be  constructed 
to  pass  through  all  three  points.  The  polynominal  over 
the  'corner'  of  Kfoif}  defined  on  the  interval  y,  ) 

can  be  and  is 

p<x)=  -IxS  |x  -  1  (90) 


Hence,  the  area  under  the  interval  from  to  X^  is 
added  onto  that  already  calculated  from  the  first  inter¬ 
val,  adding  in  the  area  (  y j)  twice.  Not  only  is 


Fig.  4  Quadratic  Interpolation  Across 
the  Discontinuity 


the  area  via  Simpson's  rule  in  error  but  note  that  by 
using  parabolas  a  larger  area  is  actually  obtained 
since  j>0<)  encompasses  the  linear  boundaries  of 
Even  for  N  odd,  parabolas  still  overestimate  the  desired 
area,  and  not  until  K0s$)is  represented  by  a  non-linear 
function  does  the  Simpson's  technique  yield  a  more 
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accurate  result  compared  to  the  trapezoid  method. 


Finite  Difference  Method 

The  method  of  finite  differences  is  a  process 
of  numerical  differentiation;  the  original  differential 
equation  is  replaced  by  finite  difference  approximations. 
This  technique  converts  the  original  differential  equa¬ 
tion  into  a  system  of  linear  algebraic  expressions 
which  can  be  solved  simultaneously  to  give  solutions 
at  a  finite  number  of  points. 

The  differential  operator  by  definition  is  ex¬ 
pressed  quantitatively  as 

,  W  y<x+Ax)-3<x) 

- Jyi -  C91  > 


However,  because  the  computer  cannot  take  the  limit 
one  must  approximate  J3.  by  using  a  small  AX  ,  that  is 


AX 


(92) 


This  quanity  is  known  as  the  first  forward  difference 
quotient.  Expression  <92  )  represents  the  change  in 
some  function  by  incrementing  X  a  positive  amount  A/. 

The  general  procedure  to  obtain  difference  ap¬ 
proximations  for  derivatives  is  to  express  several 
values  of  y  at  adjacent  points  in  a  Taylor  series 
expansion,  it  is  then  possible  to  solve  for  the  needed 
derivative  through  algebraic  manipulation. 
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To  derive  the  numerical  analogy  for  the  first 
and  second  derivatives,  expand  y  in  a  Taylor  series 
about  the  two  points  X+h  and  X-h  : 

y(x+h)=  y(x)  +  yM(*}  +  •  •  •  +  0(W*)  (93) 

y( x-h)  =  y<x)  -  V\jj'(x)  +  £ y"(x)  _  j£ y-;x)+ . ..  +  OCW1)  (94) 

Adding  equations  (93  )  and  (94  ) ,  the  sum  becomes 


y{x*h)  +  'j(x-h)  s  2  y(x)  +  \^Lj(x)  -t  •  *  •  4  O0V*)  (95) 

where  o(V)  refers  to  terms  containing  fourth  and  high¬ 
er  powers  of  h  •  By  assuming  the  higher  ordered  terms 
are  negligible  compared  to  the  lower  powers  of  V)  , 
for  small  values  of  h ,  equation  (95  )  can  be  solved 
for  yU),  and  is 

M"(x)  =  +  0(Vil)  (96) 

J  Vi 

As  V)  becomes  small  the  remainder  terms  can  be  neglected 
and  (96  )  reduces  to 


(97) 


The  above  expression  is  known  as  the  central  finite 
difference  (CFD)  approximation  to  the  second  derivative 
of  y(x)  and  has  an  associated  error  of  order  h*  . 

To  approximate  the  first  derivative,  subtract 
equation  (94  )  from  (93  )  and  obtain  the  quanity 


y  (x+h)  -  vj(x-h)  =  2h  y'<x)  +  O(h’)  (98) 
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or,  solving  for  3‘*> 

The  central  difference  approximation  for  the  first 
derivative  is  also  of  order  h*.  For  small  values  of 
h  ,vjU  becomes 


(99) 


y<x) 


&  H<*+M  ~~ 
eh 


(100) 


In  addition  to  the  central  difference  scheme  it 
is  also  possible  to  approximate  the  derivative  from 
two  adjacent  points,  a  distance h  and  2h  from  X  »  the 
technique  is  known  as  the  first  forward  difference  (FFD) 
method.  To  find  the  first  forward  difference  ex¬ 
pression  for  the  first  derivative  solve  (93  )  for 

+  o(h)  aoi) 

n 

Likewise,  for  the  second  derivative 

+  oo.)  a  02) 
h 

Note  in  the  FFD  methods  the  errors  associated  with 
equations  (101)  and  (102}  are  of  order  (■)  .  For  the  CFD 
methods  however,  the  errors  are  of  order  h*-  There¬ 
fore,  it  can  be  concluded  that  using  paints  symmetrically 
located  with  respect  to  X  give  more  accurate  results, 
by  a  factor  of  h  *  "than  those  based  on  the  FFD  tech¬ 
niques  (Ref  8«  63). 
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Matrix  Methods:  The  Integral  Form  of  Solution 


A  definite  integral  of  the  form 
written  exactly  as 

y<»=  1£.  I 


3=  J* 


can  be 


(103) 


where  the  interval  is  divided  into  N  subintervals 

of  lengths  AX,, . . . .  »<  and  XK  is  a  point  of  the  Klb 

subinterval  (Ref  Is  444).  An  approximation  to  (103 
can  be  obtained  by  expressing  ycx)  as  a  weighted  sum 
of  the  ordinates  ft**)  at  N  conveniently  chosen  points 


X,.  X2 . X*  on  the  interval  (o,k),  that  is 

£  DKf(xK)s  DNf(x„)  (104) 

K*  i 

where  is  a  weighing  coefficient  associated  with 

the  point  XK-  The  coefficient  Dr  ,  when  the  points 

*.«  . **  are  equally  spaced,  can  be  chosen  in 

accordance  with  formulas  for  numerical  integration, 
such  as  the  trapezoid  rule.  Recall  that  use  of  the 
trapezoid  rule  for  an  integral  gives 

JWixS*  ^f<Xl)4  2f<Xl)4.-  +  2f(XM.,)4f^N)^  (105) 

where  Vi=  Identifying  the  ’weighing*  coeffi¬ 

cients,  i  with  the  trapezoid  rule  it  is  found  that 

r  •  •  -i  2,2, ij  (lo6) 


r  b. 
I 


(106) 
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In  the  same  way  integral  equations  such  as 


F<*>  + 


(10?) 


can  be  approximated  in  the  form 

tJ 

yw  =  F(x>  +  Z  D*  K(x,xk);j(x><)  (108) 

where,  again  the  points  x^are  chosen  at  N  locations 
in  the  interval  (ft,fe) .  The  Devalues  correspond  to 
weighing  coefficients  based  upon  the  approximation 
method  used,  i.e.,  Simpson's  or  trapezoid  rule.  By 
requiring  that  both  sides  of  equation  (10§  be  equal 
at  each  of  the  N  chosen  points,  N  linear  equations 
results 

if 

y<*i)  =  r(xo  +  £DkK(xi,xk)^(xk)  (109) 

where  the  unknowns  y(x,) . correspond  to  the 

values  of  the  unknown  function  jjtyat  N  locations. 
Introducing  the  abbreviations 

3;=iPi),  F.xftxd,  Kij=  K(xi.x;)  (no) 

where  Kij  is  the  value  of  K(x,f)  when  X*  X;  and  §  =  Xj  , 
expression  (1C9)  becomes 


&  r  Fi  +  Kii  D«  3 1 


(111) 


The  kernel  K--  can  be  written  in  the  form  of  a  matrix 
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Also,  Dr  can  take  the  form  of  a  diagonal  matrix  by 
defining  b*  = 


quanity.  Consequently  the  integral  form  of  solution 
is  expressed,  for  computer  analysis,  as  a  problem  of 
matrix  multiplication. 


Matrix  Methods »  The  Method  of  Finite  Differences 

Just  as  was  the  case  with  the  integral  conversion, 
it  is  also  possible  to  transform  the  method  of  finite 
differences  into  a  problem  of  matrix  multiplication. 
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( 


Consider  the  general  second  order  differential  equa¬ 
tion 

'Ll  +  +  R<«)  M  H  r(X)  (H5) 

d  X*  d*  ^  ^ 

with  boundary  conditions 


y  (x,)  =  A  >  ^xn)=  B 


(116) 


Using  central  finite  differences ,  the  first  and 
second  derivatives  of  ^(x)  can  approximated  by  equa¬ 
tions  (92  )  and  (97  ) .  Substituting  these  quantities 
and  combining  terms,  equation  (115)  becomes 


4  yx+|0[i  -  r(x>V?"  (117) 

Once  again  introduce  the  notation 


E  =  >  +  ,  Rv--  ruoVf  <118) 

a 

Equation  (117  )  becomes  equivalent  to  the  system  of 
linear  equations 

JJi-iR  +  Ri  ••»*-•)  (119) 

If  the  end  conditions ,  y,*  A  and  yNr  g  »  are  substituted 
(119  )  can  be  expressed  as  the  system  of  equations 
defined  from 
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t*  A/-/ 


4  4%  4  US 

ttF3  4  S,s3  4 

4fh  4  3,®, 


w  £. .  +  w  6  +  B  E  —  Ru-i 

wn->  *m  Jam  m  h*» 


and  rewritten  in  matrix  form  as 


(120) 


Gi  Fs  0  * 

F,  G,  Ej  o 

O  F,  G,  F, 


t-i 


Rt-^5 


-  BE 


(121) 


The  left-hand  matrix  in  equation  (121)  is  known  as  a 
tridiagonal  matrix.  The  elements  on  the  principle 
diagonal,  super-diagonal,  and  sub-diagonal,  are  non¬ 
zero,  with  zero  elements  everywhere  else  (Ref  9*  104). 
Expression  (121)  can  be  written  in  the  more  concise 


matrix  form 


Bm  =  R 


(122) 


where  B  corresponds  to  the  tridiagonal  coefficient 


matrix,  and  y  and  R  are  the  perspective  matrix 
vectors.  Again,  the  differential  form  of  solution  is 
expressed,  for  computer  analysis,  as  a  problem  of 
matrix  multiplication. 


Simultaneous  Solution  of  Linear  Equations  ^ 

The  solution  for  a  system  of  equations,  A  x  =  \> 
is  most  efficiently  solved  by  eliminating  the  unknowns, 
and  the  most  commonly  employed  method  is  the  Gauss 
elimination  process.  Consider  as  an  example  the  follow¬ 
ing  set  of  four  equations! 


c 

c 


II 

3‘ 

Cja. 

C„ 

C-JH 

HI  . 

<^H3 

C-HH 

n 

*3 

-1 

(123) 


The  solution  vector,  X  ,  is  the  desired  quanity  and 
will  be  obtained  through  algebraic  manipulation  of 
the  coefficient  matrix  in  (123).  During  the  process 
of  this  matrix  manipulation  it  should  be  remembered 
that  the  solution  of  %  remains  unchanged  if  any  of 
the  following  operations  are  preformed: 

(1)  Multiplication  or  division  of  any  equation 
by  a  constant. 

(2)  Replacement  of  any  equation  by  the  sum  or 
difference  of  that  equation  and  any  other 
equation. 

Gauss  elimination  is  simply  a  sequential  app¬ 
lication  of  the  row  operations  (1)  and  (2)  above.  In 
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the  algorithm,  the  top  equation  is  first  divided  by  C,, 


(124) 

where  the  primes  denote  elements  whose  values  have 
been  changed.  The  first  equation  is  now  multiplied  by 
and  subtracted  from  the  second  to  eliminate  element 
C-21  .  Likewise  the  first  equation  can  also  be  multi¬ 
plied  by  C31  and  subtracted  from  the  third  equation, 
then  multiplied  by  Cm  and  subtracted  from  the  fourth. 

In  this  manner  the  entire  first  column  below  Cu  has 
been  cleared  to  zero  and  the  set  appears  as 


) 

<4 

<» 

ST* 

cti 

Co 

cn 

C,x 

Cj«*  *3 

r* 

Cm 

CHJ 

Ji 

- - 

—  — 

/ 

1 

c; 

r  '  c  ' 

>-13  '“14 

*  \ 

t 

0 

Cn 

/  1 
Cjj  Cm 

*z 

rt 

0 

ci 

c»  ei. 

* 3 

1 

r3 

0 

cw  CJ<1 

XI 

1 

(125) 


During  the  column  manipulation,  the  first  row  is 
termed  the  pivot  row  and  C„  the  pivot  element.  The 
second  row  now  becomes  the  pivot  row  and  C ^  "the  pivot 
element.  The  second  equation  is  divided  by  to 
make  the  main  diagonal  element  1.  Multiplication  of 
the  second  equation  by  cn  and  subtraction  from  the 
third,  and  then  multiplication  by  and  subtraction 
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from  the  fourth,  clears  the  second  column  below  the 
main  diagonal  to  zero.  Similar  operations  with  the 
third  and  fourth  rows  as  pivot  rows  finally  yield 


1 

c* 

c  # 

'■'13 

r * 

0 

\ 

c  * 

0 

O 

\ 

C3M 

0 

0 

0 

(126) 


where  the  stars  indicate  elements  which  have  been 
modified  several  times  from  their  original  values. 
Gauss  elimination  is  sometimes  called  triangular- 
ization  because  the  coefficient  matrix  is  upper  tri¬ 
angular,  i.e.,  elements  above  the  main  diagonal  are 
non-zero  while  those  below  equal  zero.  The  bottom 
equation  in  (126)  now  directly  gives  the  value  of 


The  third  equation  is 

x3  *e  X4  =  *3 

Because  X4  is  known,  Xj  can  be  solved  to  yield 


(127) 


(128) 


(129) 


Repeated  back  substitutions  will  evaluate  one  new 
unknown  from  each  new  equation.  The  unknown  vector  x 
will  have  been  completely  determined  when  the  top  equa- 
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tion  is  solved  for  X(. 

It  should  be  mentioned  that  one  computational 
difficulty  can  arise  with  the  Gauss  elimination  tech¬ 
nique.  The  pivot  element  in  each  row  is  the  element 
on  the  main  diagonal.  By  the  time  any  given  row  be¬ 
comes  the  pivot  row,  the  main  diagonal  element  in  that 
row,  i.e.,  the  pivot  element,  will  have  been  modified 
from  its  original  value  several  times  over.  Under 
certain  circumstances,  the  diagonal  element  can  be¬ 
come  very  small  in  magnitude  compared  to  the  rest  of 
the  elements  in  that  row,  as  well  as  inaccurate. 

This  can  result  in  an  erroneous  solution  vector.  The 
problem  can  be  effectively  treated  by  interchanging 
the  columns  of  the  matrix  to  shift  the  largest  ele¬ 
ments  (in  magnitude)  in  the  pivot  row  into  the  diag¬ 
onal  position.  The  largest  element  then  becomes  the 
pivot  element.  With  each  new  pivot  row  the  operation 
is  repeated.  This  scheme  is  known  as  partial  pivot¬ 
ing  and  minimizes  the  roundoff  errors. 

In  addition,  it  is  important  to  note  the  two 
different  types  of  coefficient  matrices  which  can 
be  encountered.  For  the  integral  solution  the  system 
of  N  algebraic  equations  will  have  a  'dense'  NxN  co¬ 
efficient  matrix,  few  zero  off  diagonal  elements. 

The  solution  of  algebraic  equations  for  the 
finite  difference  technique  yields  a  banded  coefficient 
matrix  and  is  'sparse',  many  zero  off  diagonal  ele- 
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ments.  These  sets  commonly  arise  in  the  numerical 
solution  of  partial  differential  equations  and  in  the 
use  of  finite  differences.  Solutions  of  banded  systems 
are  considerably  easier  than  those  for  full,  dense 
coefficient  matrices. 

As  an  example,  consider  this  set  having  a  tri¬ 
diagonal  coefficient  matrix: 


tj  c,  o  «  •  • 

x, 

~ 

a*  bt  c*  0  •  • 

r* 

O  a5  b»  c3  0 

X3 

, 

•  t 

• 

•  • 

• 

• 

0 

*  k-l  Sl-I 

XN., 

'm-\ 

_r»L 

(130) 


where  the  main  diagonal  elements  are  denoted  as  b  ,  and 
the  diagonal  elements  below  and  above  the  main  diagonal 
are  Q  and  G  respectively. 

Applying  Gauss  elimination  to  (130  )  only  one 
element  (one  of  the  Cl's)  will  be  eliminated  in  each 
column  because  all  remaining  entries  below  the  main 
diagonal  are  zero.  Also,  no  entries  outside  the  tri¬ 
diagonal  band  are  changed  from  zero  in  the  course  of 
the  elimination  process.  After  the  bottom  row  has 
been  reached  (130  )  becomes 
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coefficient  matrix  since  only  the  three  diagonal 
vectors  Q.  ,  b,  and  C  are  required  (Ref  10«  90-98). 
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III.  .One-Dimensional  Problems 


Computer  Programs 

The  computer  programs  written  for  all  one-dimen¬ 
sional  cases  were  very  similar;  the  programs  used 
common  expressions,  symbols,  and  followed  parallel 
sequences  so  that  minimum  alterations  were  necessary 
to  adapt  each  routine  to  new  cases.  The  flow  diagrams 
for  both  numerical  integration  and  differentiation 
appears  in  Figure  5 . 


Read  N,  number  of  data  points 
BC1,  boundary  condition  1 
BC2,  boundary  condition  2 _ 


Fill  X  array,  X,»  X? . *  Xm 

Generate  ^(x  4)  or 

Fill  _ _ 


Perform  matrix  multiplication  and  all 
other  operations  to  generate  the 
coefficient  matrix 


CALL  LEQT1F  Subroutine 

International  Mathematics  Science  Library 
(IMSL)  -  Solves  a  linear  system  of  equations 


Print  out  solution! 
Ivector 


Fig.  <5  Flow  Diagram  for  Computer  Programs 

Using  Numerical  Integration/Differentiation 
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Because  the  system  of  matrix  equations  are  very 
similar,  i.e.,  expressions  (114)  and  (122),  the  only 
modification  between  the  two  programs  were  in  the 
coefficient  matrices.  For  numerical  integration,  the 
coefficient  array  was  computed  from  matrix  multi¬ 
plication  with  the  Green's  kernel,  and  the 

weighing  coefficient,  D^*  In  the  finite  difference 
method  however,  the  coefficient  array  was  filled 
directly  because  of  the  simple  band  structure  simply 
by  noting  the  sequential  arrangement  of  internal 
terms . 

LEQT1F  -  IMSL  Subroutine 

The  International  Mathematics  Science  Library 
(IMSL)  routine  LEQT1F  was  used  to  solve  all  linear 
algebraic  equations  generated  in  the  matrix  form¬ 
ulations.  LEQT1F  solves  the  set  of  linear  equations 
fix  =  B  for  X  ,  given  the  NxN  matrix  and  Nxl 
matrix  B  .  The  solution  X  will  be  the  exact  solution 
without  any  roundoff  error.  If  such  a  solution  can¬ 
not  be  obtained  a  warning  is  given.  LEQT1F  performs 
Gauss  elimination  and  uses  partial  pivoting  of  the 
array  elements. 

Case  I 

The  first  problem  addressed  in  this  thesis,  ex¬ 
pressed  mathematically  is 


to  =  - 

<Jx* 


3Cx) 


(132) 


where  y<x>  is  the  inhomogeneous  term,  with  boundary- 
conditions 


y(o)  = 


=  O 


=  i 


(133) 


Analytical  Solution 

Equation  (132  )  is  equivalent  to 


+  b‘*>  -  ° 


ci** 

and  can  be  solved  in  a  straightforward  manner  by 
assuming  yOO  to  be  a  linear  combination  of  the  ex¬ 
ponential  function,  that  is 

ij(x)=  eDx 


(134) 


(135) 


Dx 


Substituting  (135  )  into  (134 )  and  dividing  by  C 
yields 

\  =  O  J  Ds  ±  4  (136) 

Using  the  two  roots  for  0  ,  the  solution  to  (134 )  be¬ 
comes 

-  fte+,Xf  Be  11  =  /{cos  x  ■>  B'smx  (137) 


vjOO 


By  applying  the  boundary  conditions,  y(O)r-0and  1)01=1, 
the  constants  ft*  and  g/  are  evaluated.  The  complete 


solution  is  found  to  be 


smx 
sin  | 


and  is  graphed  in  Figure  6. 


(138) 


Integral  Solution 

To  convert  differential  equation  (132  )  into 
integral  form,  first  integrate  both  sides  of  (132) 
with  respect  to  X  over  the  interval  (QX) 


(139) 


or 


#<x>  =  C  -  fy-oa* 


(140) 


where  the  constant  C  is  from  the  lower  limit, 

A  second  integration  over  the  same  interval  and  using 
equation  (6)  for  the  right-hand  side,  equation  (140  ) 
becomes 

* 

(141) 


y(X)  =  CX  - 


The  constant  C  can  be  evaluated  from  the  boundary 

condition  Jj(0*|and  equation  0-41  )  takes  the  form 

I 

x  +  J"  K(x,f)  <142) 

where  the  kernel  is  defined  as  the  function 
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3(0)  =0,  «j(l)  =  I 

Hin)5  Sinx 


Analytic  Solution,  Case 


f  *<»-§> 

[s(>  -X) 


x<% 

x>S 


(143) 


Rearranging  terms,  (132  )  becomes  equivalent  to  the 
expression 

r« 

M<x)  -  J  =  x  U44) 

Now,  use  the  approximation  for  trapezoid  integration 
for  the  integral  and  (144  )  becomes 

y(x)  -  2.  0^K(x,xk)^xk)  =  x  (145) 

Expressing  and  as  matrices,  equation  (145  )  is 

transformed  into  the  system  of  equations 

-  u  =  x  u«) 

or 

(i -  K  t>)  jj  =  x  <14?> 

For  the  case  where  the  interval  (0,1)  is  divided  into 
four  equal  partitions,  N=5»  at 

X,  -  O,  X3*'/2;  =  Xs=  1 

the  matrix  equation  0-47  )  takes  the  form 


3 
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ih~  ° 

ft  it  vh 

y,-k%=  '/i 

CH  ~  ife  ^  fe* i  Jfo*  ^ 

V  1 

Solution  of  this  set,  rounded  to  six  decimals,  is 

\Ji  =  °»  ^i=  -294274,  cj3=  .570156,  t/4=  .810403,  y5=  1 

The  true  values  obtained  from  the  analytical  solution, 
equation  (13  8),  are 

%  s  °»  Jjk*  .294014,  ^3=  .569747,  .810056,  1 


Finite  Differences 

In  a  similar  manner  as  was  done  to  the  integral 
equation,  the  method  of  finite  differences  is  now 
applied  to  the  original  differential  equation.  By 
substituting  the  second  derivative  CFD  quotient,  equa¬ 
tion  (97),  back  into  the  original  differential  equation 
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forgot),  and  collecting  terms,  equation  (132)  takes 
the  form 

&-»  4  +  &♦»*  °  (150) 

where  N  is  the  number  of  points  over  the  interval  (0,1) 
and  h ,  the  interval  spacing.  A  set  of  linear  equa¬ 
tions  can  be  generated  from  (150  )  by  letting  L  run 
through  its  sequence  up  to  the  value  N-l.  For  N=5  and 
using  the  same  five  points  as  before,  the  central 
finite  difference  method  will  give  the  following  set 
of  linear  equations! 


Hi 

'8$ 

=  O 

% 

■8  a, 

+  S4 

=  O 

&  ~ 

=  O 

in  matrix 

form 

1 

-2L 

lb 

\ 

0  0 

— i 

0 

0 

1 

O 

-31 

lb 

1  “ 

1  0 

2 L  1 
lb _ _ 

1 - 

r 

1  O  O 

(15D 


(152) 


Using  the  initial  boundary  conditions,  ^rOand  |  , 
the  following  values,  correct  to  six  decimal  places, 
are  obtained  for  y*)' 

^1=  °»  Jt=  -294274,  .570156,  .810403,  %=  1 
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Note  that  the  same  values  were  obtained  using  the 
trapezoid  rule  for  numerical  integration. 


First  Forward  Difference 

Besides  the  method  of  central  finite  differences, 
the  first  forward  difference  expression,  equation  (1025, 
was  substituted  into  (132  )  for  a  comparison  of  the 
two  difference  methods.  Equation  (132  )  now  becomes 

y-(l+h)  t  $+2  =  o  L-  ».*, . N‘2‘  U53) 

Again,  for  the  same  five  points,  (153  )  generates  the 
linear  equations 

U'jrt'k+h  s° 

Slfe'2***,  =°  (l54> 


or  the  equivalent  banded  matrix 


(155) 


This  set  of  equations,  correct  to  six  decimal  places, 


has  the  solution 


y,  =  o,  yt=  .266667,  %=  .533333,  .783333,  J5=  1 

For  a  comparison  of  the  relative  errors  associated  with 
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each  of  the  two  difference  methods  see  Figure  7*  The 
relative  error  data  is  displayed  in  Table  I. 


with  the  particular  method  being  employed.  All  other 
error  sources,  such  as  truncation  and  internal  round¬ 
off,  due  to  the  arithmetic  process,  are  insignificant 
when  compared  to  the  errors  created  by  a  specific 
routine.  By  neglecting  these  other  errors  it  will  be 
possible  to  predict  how  large  an  error  should  occur  and 
compare  this  to  that  actually  observed. 


Integral  Equations 

The  error  bound  using  the  trapezoid  rule  (Ref  12 i 
405-408)  is  given  by  the  expression 

.  I  >  .  .  t  9  I  A  ^  _  I 


-  T 


(156) 
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X  *80883  3AIi«“l38 


where  ^f(Oci§  is  the  required  integral 

T  -  number  calculated  via  trapezoid  rule 

(b-a)  _ 

the  boundary  interval 
Vl  -  increment  size 

fA  -  upper  bound  satisfying  the  condition 


A\  >  |f?o 


(157) 


To  predict  how  accurate  the  trapezoid  rule  is, 
the  upper  bound  M  must  first  be  found.  Again,  it  will 
be  assumed  that  the  only  error  in  the  integral  solu¬ 
tion  of  equation  (144  )  is  that  due  to  numerically 


approximating  the  integral: 

i 


(158) 


which,  inserting  the  appropriate  limits,  becomes  equiv¬ 


alent  to  two  integrals: 


C  c‘ 


(159) 


The  error  will  be  computed  at  the  point  X=)£  since  a 
number  of  values  from  Table  I  exist  for  comparison. 

Substituting  X*  expression  (159)  becomes 

*/2  1 

jsja’is  + 

The  function  ,  from  expression  (156),  can  be 
taken  as  the  sum  of  both  integrands  above;  that  is 


(160) 


■fa)*  ♦  i 


(161) 
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where 


f,(o=  *«<*  f;<o=  fl-pyo  (162) 

^  2. 

Taking  the  second  derivative  of  f<5>  and  substituting 
this  quanity  back  into  (157  )  will  yield  the  bound  M . 

It  is  found  that 

f  (§)  = 

The  bound  A\  is  therefore 

M  >  l  ■(  "(5)1  = 


(163) 


(164) 


The  function  ys)  has  its  maximum  value  at§=l  and  y(f) 
has  the  value  of  1,  from  Figure  6.  Therefore,  A1 
will  be  the  upper  bound 

*  Yz 


The  total  error  at  the  point  using  trapezoid  rule, 

becomes 


Total  error  = 


(b-a)MW1 

\l 


(166) 


As  h  is  decreased  the  associated  error  should  also 
decrease  by  a  factor  of  H*.  This  trend  is  indeed  ob¬ 
served  in  the  integral  solution  values  shown  in 
Table  II. 


Simpson’s  Comparison 

In  addition  to  use  of  the  trapezoid  rule,  Simpson's 
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Table  II.  Error  Trends  Using  Trapezoid  Rule,  x=.5 


Increment  Size 
h 

y(x) 

Exact 

y(x) 

Calculated 

Actual 

Error 

Calculated 

Error 

1/2 

.569747 

.571429 

-.001682 

.0104167 

1/4 

.569747 

.570156 

-.000409 

.0026042 

1/6 

.569747 

.569928 

-.000181 

.0011574 

1/8 

.569747 

.569849 

-.000102 

.0006510 

1/50 

.569747 

.569750 

-.000003 

.0000167 

rule  for  numerical  integration  was  also  used  on  (134). 

An  interesting  though  not  unpredictable  trend  was  noticed 
in  the  data,  Figure  8. 

For  %Ot  the  computed  array  elements  corres¬ 
ponding  to  the  even  intervals,  i.e.,  odd  array  numbers 
were  Simpson's  method  can  be  used,  showed  much  smaller 
error  values  as  compared  to  the  even  array  elements. 

Also,  as  the  array  elements  increased  the  relative 
error  values  decreased  over  both  odd  and  even  elements. 
The  relative  error  is  high  at  the  even  array  elements 
due  to  the  additional  sum  from  an  extra  partition  near 
the  'corner'  in  the  quadratic  interpolation  scheme. 

The  effect  is  to  change  the  true  area  by  an  amount  of 
the  previous  partitions  area.  The  relative  error 
values  decrease  because  the  partition  areas  decrease 
after  the  middle  value ,  For  the  lower  points  on 

Figure  8,  where  Simpson's  rule  is  legal,  the  parabolas 
tend  to  approximate  the  kernel  closer  and  have  less 
overlap  across  the  'comer'  sections  of  KM) .  The 
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SOPTS.  SlttSO*  RUL£ 


Using  Simpson  • 


SOFTS.  TRRPEZOIO  RLE 


Fig.  9  Relative  Errors  Using  Trapezoid  Rule 


relative  errors  for  the  same  array  elements  using 
the  trapezoid  rule  is  shown  in  Figure  9. 

Finite  Differences 

The  major  errors  associated  with  the  finite 
difference  methods  are  from  truncations  in  the  Taylor 
series  expansions  used  for  y(x).  These  errors  can  be 
reduced  by  taking  more  terms  in  the  series  expansion, 
however  doing  so  does  not  necessarily  guarantee 
numerical  stability.  Because  of  these  facts  it  is 
extremely  difficult  to  predict  the  relative  size  of 
error  which  should  result  from  using  the  finite  differ¬ 
ence  method.  The  error  is  carried  in  all  linear 
equations  generated  by  the  recurrsion  expression;  when 
the  number  of  equations  becomes  large,  the  resulting 
error  combinations  become  exceedingly  complicated. 

It  is  possible  however  to  note  the  error  trends  in 
both  the  CFD  and  FFD  methods. 

Recall  that  the  CFD  method  was  of  order  Mb').  terms 
of  order  H*  and  higher  were  neglected  in  the  central 
difference  quotient,  the  second  derivative  became 


It  is  reasonable  to  assume  that  errors  associated  with 
this  method  are  as  a  direct  result  from  neglecting 
the  truncated  terms  of  order  h*  and  higher;  there¬ 
fore,  the  trend  of  error  should  be  of  order 


The 


error  trends  using  the  central  and  first  forward 
difference  quotient  are  compared  in  Table  III. 


Table  III.  Error  Trends  Using  Finite  Difference 
Methods  for  y"+y=0,  x=.5 


Increment  Size 
h 

y(x) 

Exact 

y(x) 

CFD 

Error 

y(x) 

FFD 

Error 

1/2 

.569747 

.571429 

-.001682 

.500000 

.069747 

1/4 

.569747 

.570156 

-.000409 

.533333 

.036414 

1/6 

.569747 

.569928 

-.000181 

.545455 

.024292 

1/8 

.569747 

.569849 

-.000102 

.551576 

.018171 

1/50 

.569747 

.569750 

-.000003 

— 

Note  from  Table  III,  that  as  the  interval  spacing  is 
halved,  from  to  ln« ^ ,  the  error  for  the  CFD  method 
decreases  on  the  order  of  or  as  In1.  For  the  FFD 
method,  the  error  decreases  on  the  order  of  h  or 
approximately  by  Vfc  . 

Times  of  Solution 

The  method  of  central  finite  differences  was 
compared  to  the  integral  equation  method  for  speed. 
The  data  from  this  analysis  is  contained  in  Table  IV. 

Table  IV.  Computer  Time  -  y"+y=0 


Method  Compilation  Time  Execution  Time  Usage 
_ ( sec) _ (sec) _ (Kilo-word) 

Integral  0.735  1.178  199.410 

Equation 

CFD  0.531  0.5^6  166.875 
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The  results  of  this  comparison  revealed  that  the  CFD 
method  required  53*6$  less  execution  time,  was  com¬ 
piled  27.8$  faster,  and  used  less  storage  than  the 
integral  equation  method.  The  time  savings  can  he 
attributed  to  the  different  coefficient  matrices  gen¬ 
erated  by  each  method.  The  CFD  method  gave  a  tridiag¬ 
onal,  sparse,  49  x  49  coefficient  matrix.  The  integral 
method  however  resulted  in  a  full,  50  x  50  coefficient 
matrix,  which  required  more  time  and  more  manipulations 
to  solve  for  the  solution  vector  MIX)  . 


Accuracy 


The  central  finite  difference  method  was  com¬ 


pared  to  the  integral  method  for  accuracy.  Both 
methods  were  found  to  be  equivalent?  that  is,  the  ex¬ 
act  same  solutions  were  obtained  from  both  procedures, 
identical  to  12  decimal  places.  Equivalency  of  the  two 
methods  was  verified  through  matrix  multiplication  of 
equations  (  142)  and  (150).  This  was  accomplished  by 
rewriting,  in  corresponding  matrix  form,  equation 
(147)  and  (152).  The  integral  equation  was  expressed 
symbolically  as 


/lx=  F 


(168) 


where  n  represents  the  coefficient  matrix,  X  the 
solution  vector,  and  F  some  function  vector.  The 
unknown  X  ,  can  be  solved  and  is 


62 


(169) 


X  8 


where  n  is  the  matrix  inverse  of  A  .  In  a  similar 
manner  the  CFD  method,  equation  (150),  was  expressed 
symbolically  as 


=  R  (170) 

where  b  represents  the  coefficient  matrix  associated 
with  (121).  The  unknown  y  ,  can  be  solved  by  matrix 


inversion  and  found  to  be 

T'l 


(171) 


Since  both  methods  yield  identical  solutions,  that  is 


x  *  y  (172) 

the  expressions  on  the  right  side  of  equations  ( 169) 
and  (171)  should  also  be  equivalent,  therefore 

ft  'f  1  ft  R  <1?3> 

The  necessary  condition  for  equivalency  is 

f  -  tm) 

Equation  ( 174)  is  easy  to  verify  since  F  and  A  are 

m  sr 

already  known  and  the  coefficient  matrices,  Q  and  B  * 
are  internally  generated.  Using  7  iterations,  equa- 


< 

4 


tion  (17*0  was  verified  by  substituting  the  known 

y  s./  —  — 

quantities  n  ,  B  ,  and  R  ,  the  calculated  F  was 
found  to  be 

(175) 

The  true  value  of  F  ,  from  the  integral  equation 
computer  program,  was 

F  = 

The  reason  for  the  discrepancy  in  the  sixth  decimal 
place  pertains  to  a  rounding  error  in  the  multipli¬ 
cation.  The  original  computer  program  uses  all  14 
decimal  places  for  internal  computations.  In  this 
example  only  six  significant  figures  were  used  and  a 
rounding  error  results  in  the  sixth  decimal. 

Conclusions 

It  was  unforseen  that  both  the  CFD  and  integral 
approach  would  be  identical.  Though  the  two  methods 
were  solved  using  the  methods  of  linear  algebra  each 
was  characterized  by  a  distinct  different  coefficient 
matrix,  one  being  sparse  for  the  CFD  method,  and  the 
other  being  full  for  the  integral  approach.  Both 
methods  were  shown  to  be  of  the  same  order  of  error, 


.142857 

.285714 

.428571 

.571429 

.714286 

.857143 


(176) 


.142857 

.285713 

.428571 

.571428 

.714285 

.857143 
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h  .  It  is  not  possible  to  generalize  at  this  point 
any  other  conclusions  except  that  both  methods  are 
equivalent.  The  CFD  solutions  were  noted  to  be  com¬ 
puted  faster  in  all  circumstances  but  only  because  of 
the  tridiagonal  coefficient  matrix.  All  computer  solu¬ 
tions  proved  to  be  accurate  with  no  significant  round¬ 
off  errors. 


Case  II 

The  second  problem  investigated  was  the  one-dimen¬ 
sional  inhomogeneous  Helmholtz  equation  with  homogeneous 
Dirichlet  boundary  conditions.  Expressed  mathematically, 
the  problem  is 

+  ^j<*>  =  i  u??) 


where 


,  with  boundary  conditions 


u(o)=  O, 


(178) 


Analytical  Solution 

As  with  all  inhomogeneous  2nd  order  differential 
equations  the  solution  can  be  expressed  as  a  linear 
combination  of  solutions:  homogeneous  +  particular. 

The  complete  general  solution  to  (177)  is 

y(x)s  flcos  Ax  +  B  Sin  Ax  +  ^  (179) 

Using  boundary  conditions  (178),  the  constants  A  and  B 
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are  evaluated,  the  unique  solution  is 

I?  x  cos  |p  *  4  1 ) 


and  is  graphed  in  Figure  10. 


(180) 


Integral  Solution 

The  differential  equation  can  be  converted  to 
integral  form  by  first  rearranging  terms  and  solving 
for  the  second  derivative.  Equation  (177)  becomes 


is.  -  i  -  x1 


(181) 


Integrate  each  side  over  the  interval  (o,x)  to  obtain 


= c  +Jo-/y*>)<Jx  (1 

*  o 

where  the  constant  C  corresponds  to  .  Integrate 

a  second  time  over  the  same  interval,  equation  ( 182) 


(182) 


becomes 


rY 

y(x)»  cx+J 

o  o 

Using  (6)  on  the  double  integration  gives 

X 

yfx)=  Cx  -t 


(183) 


(184) 


The  constant  C  can  now  be  evaluated  by  using  the  second 


boundary  condition,  ry- 0 ' 

i 

jo)- o  -  c  + 


(185) 
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Fig.  10  Analytic  Solution,  Case  II 


(186) 


Solving  for  C 

i 

c  =  Jo-DOtyo-Ol* 

o 

Substituting  the  expression  for  C  into  equation  ( 184) 
yields 

I 

yx)=  <4  (187) 

o 

with  the  kernel  defined  by  the  expression 


y<  $ 

Sh-x) 

v>§ 

(188) 


Equation  (18?)  represents  the  equivalent  integral  form 
for  the  one-dimensional  Helmholtz  equation. 

Equation  ( 1 87 )  can  be  rearranged  by  multiplying 
through  by  K<«)  and  separating  the  integral.  The  in¬ 
tegral  equation  containing  y  becomes 
i  1 


y*)  -JVfctf)  AyoJ*  =  -f  Kfr,§)i$ 


(189) 


The  integral  on  the  right  simplifies  by  replacing 
with  its  defined  values  over  the  appropriate  limits 
of  integration,  that  is 

o 


(190) 
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Integrating  out  the  S  variable  on  the  right  side, 

equation  (190)  becomes 

l 

yc*)  -  JVx,g)A*ys)cls  =  f(x-0  (191) 

o 

The  trapezoid  approximation  can  now  be  used  for  the 
integral  and  equation  ( 1 9 1 )  becomes  equivalent  to 

y(x)  -  2_  =  X-(x-i)  (192) 

Expressing  Dj^  and  as  matrices,  equation  (192) 

transforms  into  the  system  of  linear  equations 

(i  -  -  x  Cx-0  (iM) 

and  is  solved  by  the  method  of  Gaussian  elimination. 


Green's  Integral 

In  addition  to  the  integral  representation  of 
equation  ( 177)  it  is  also  possible  to  express  the 
solution  of  jfM  by  the  Green’s  integral,  equation  (27) 
From  Appendix  B  ,  the  Green's  function  for  the  Helm¬ 


holtz  operator  was  found  to  be 


G 


^sm  sm  A  x 

x  <  f 

A  sm  XSL 

sin  A (t-x)  sin  X$ 

x  >% 

Asm  AX 

(194) 


Therefore  the  equivalent  expression  for  the  Green's 
integral  solution  for  is 
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(195) 


y(x)  -  -  j"  G(x,§)^S 


Using  the  trapezoid  rule  equation  ( 195)  becomes 

yx>--  - 1  \zu,yk) 


(196) 


K- I 


or,  expressed  in  matrix  notation 


(197) 


where  D  is  the  column  vector  of  dimension  Nxl  and 
G  the  Green's  matrix,  defined  by  expression  (194), 
of  dimension  NxN 


Finite  Differences 

The  method  of  CFD  was  used  to  numerically 
differentiate  equation  (177).  Substituting  the  central 
difference  quotient,  equation  ( 97  ) ,  for  the  second 
derivative  and  collecting  common  terms,  equation  ( 177) 
becomes 

tj.  -(2 s  h  ;* 2.3,-. •, N~*  (198) 

The  equivalent  matrix  expression  has  the  tridiagonal 
form 
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I  I  o 

O  l  -(2-Xlh*)  I 


i  -  fr-A'h*)  l  O 
0  I  | 


(199) 


( 


Using  the  boundary  conditions,  yf=  O  and  yM=  O  , 
expression  (  199)  can  be  written  in  the  form 


(2 -XT)  i  o 
i  -a- xt)  i 


i  -a-w)  i 
>  -a-yif) 


f 


If 


(200) 
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The  dimensions  of  all  three  matrices  have  been  re¬ 
duced  to  N-2  .  In  more  concise  form  equation  ( 200) 
becomes 


(201) 


where  $  corresponds  to  the  band  coefficient  matrix 
and  R  the  function  vector  defined  by  the  constant 


Relative  Error 

For  a  comparison  of  the  relative  errors  asso¬ 
ciated  with  the  three  numerical  methods  based  on  the 
Helmholtz  equation,  see  Figure  11.  The  relative  error 
of  the  three  numerical  methods  is  shown  in  Table  V. 
Once  again  identical  results  were  obtained  using  the 
integral  equation  and  CFD  method 


Table  V.  Relative  Error  {%)  -  Computer  Solutions 
for  y"+)fy=l,  x=.5 


Number  of  Iterations 
_ per  Interval 


Integral-CFD  Green's  Integral 

Method  Method 


2 

4 

6 

8 

50 

Error  Analysis 


35.253527 

7.609482 

3.472912 

1.971673 

0.051070 


51. 201617 
11.842646 
5.194057 
2.908348 
0.074035 


The  same  assumptions  will  be  used  for  the  error 


54.47*]  t*  t  K  si  i  9^4 

j  0  INTEGRAL  -  CFG  HETH. 

J  ?  Q  GREEN'S  FUN/OIFF-OT. 

1  /  CONPRR I SON/ERROR 
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analysis  as  were  presumed  in  the  first  case  study, 
all  errors  are  strictly  due  to  the  method  used  and 
errors  generated  in  the  computational  process  could 
be  neglected. 


Integral  Equation  Errors 

The  integral  equation,  equation  ( 1 87 )  and  the 
Green's  integral  equation,  equation  (195),  have  both 
been  shown  to  be  equivalent;  therefore,  they  should 
have  associated  with  them  the  same  order  of  error 
because  the  trapezoid  rule  is  used  to  evaluate  each 
integral.  To  predict  the  order  of  error  the  upper 
bound  must  first  be  found.  The  Green's  integral 
equation  will  be  used  for  this  analysis.  Equation 
(195)  can  be  written 

fefr  *  x  <202> 

o  X 

where  the  expression  for  GCx,  f)  has  been  used  over  its 
defined  limits.  The  error  will  again  be  computed  at 
the  point  X*  ^2  ,  using  the  values  and  X  -  1 


Equation  (202)  becomes 


ijoo=  * Ai  4  J 5,0 ¥<'-«> j 

The  total  error  of  using  expression  (203)  will  be  no 
greater  than  the  sum  of  errors  from  each  integral 

-.is£|s.n?n|  + 


(203) 


(204) 
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Taking  the  first  two  derivatives  of  (204)  and  using 
the  expression  for  the  error  bound,  M  is  found  to  be 


A1  ^  VT h>l  r  +^Cos3f%J 


*  'J 


(205) 


The  maximum  value  for  ^Svn^jS^  cos^JTfJ  occurs  when 


5  = 


-  jr 


,  therefore  the  bound  for  M  is 


M=  (.\5)^5l/’|.37S3)  =  V.r?// 

And  the  associated  error  is  of  the  order 


(206) 


£  t't-or  ^  (V»  — <».) Z'AVi*"  =  ,3%2L  b2-  (207) 

12. 

As  h  is  decreased  the  associated  error  also  decreases 
as  a  factor  of  .  This  is  observed  in  Table  VI. 


Table  VI.  Error  Trends  -  Computer  Solutions 
for  y"+Axy=l,  x=.5 


Increment  y(x) 
Size,  h  Exact 

y(x)  Cal 
Int-CFD 

Error 

y(x)  Cal 
Green's 

Error 

1/2 

.108716 

.070390 

.038326 

.053052 

.055664 

1/4 

.108716 

.100443 

.008273 

.095841 

.012875 

1/6 

.108716 

.104940 

.003776 

.103069 

.005647 

1/8 

.108716 

.106572 

.002144 

.105554 

.003162 

1/50 

.108716 

.108660 

.000056 

.108635 

.000081 

In  addition  to  the  above  table  the  relative  error  was 
examined  at  each  array  element  when  the  interval  (0,1) 
was  subdivided  into  50  partitions.  Because  the  integral 


and  CFD  methods  gave  identical  results  their  errors 
both  appear  in  Figure  12.  The  Green's  integral  was 
also  analyzed  at  each  array  element  and  appears  in 
Figure  13.  It  can  be  seen  from  both  sets  of  dats  that 
the  errors  are  symmetric  and  largest  in  the  neighbor¬ 
hoods  of  the  end  points.  This  is  attributed  to  the 
functional  form  of  both  3<x>  and  the  Green's  function. 
At  the  end  points  the  two  functions  have  their  steep¬ 
est  incline  therefore  larger  errors  are  introduced 
over  the  same  size  increments  across  these  regions. 

At  those  points  along  the  center  of  the  interval, 
where  the  functions  change  more  slowly,  less  error  in 
the  increment  is  introduced.  Note  that  the  Green’s 
integral  method  is  more  accurate  at  the  outer  elements 
(1)  -  (10)  and  (40)  -  (49).  The  integral  -  CFD  method 
is  more  accurate  at  the  inner  elements  (11)  -  (39). 

Also  the  error  trend  in  the  Green’s  method  is  averaged 
out  to  almost  a  constant  . 0 75%  while  the  CFD  -  integral 
method  appears  parabolic  across  the  region. 


Finite  Differences  Errors 

No  additional  error  analysis  was  carried  out  for 
the  method  of  finite  differences.  It  was  noted  that 


errors  associated  with  the  CFD  method  once  again  were 
of  order  f)  due  to  the  second  order  truncation  in  the 


Taylor  series  expansion  for 


Timpfi  nf  Solution 

The  two  numerical  untegration  routines  were  com¬ 
pared  to  the  CFD  method  for  speed  and  storage.  The 
interval  (0,1)  was  divided  into  50  partitions,  re¬ 
sulting  in  a  set  of  49  linear  equations  for  the  CFD 
method  and  a  set  of  51  linear  equations  for  the  inte¬ 
gral  method.  Using  the  Green's  integral  expression 
the  solution  was  obtained  by  matrix  multiplication. 
The  data  for  these  three  methods  are  contained  in 
Table  VII. 

Table  VII.  Computer  Time  -  y"+  \y=l 


Method 

Compilation  Time 
( sec) 

Execution  Time 
(sec) 

Usage 

(Kilo-word) 

Integral 

Equation 

.784 

1.177 

194.460 

Green’s 

Integral 

.600 

0.441 

116.319 

CFD 

.548 

0.558 

170.556 

The  CFD  method  was  2.1  times  faster  than  the  integral 
equation  method.  In  addition,  it  was  compiled  in  30.1% 
less  time  and  used  12.3%  less  central  storage.  The 
Green's  integral  method  was  found  to  be  quickest  and 
used  the  least  amount  of  storage. 

Conclusions 

The  CFD  and  integral  equation  methods  both  proved 
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i 

t 


( 


to  be  identical.  Though  the  two  were  again  charcter- 
ized  by  distinct  coefficient  matrices,  they  were  found 
to  be  equivalent.  A  second  integral  approach  was  ana¬ 
lyzed  by  solving  for  the  exact  Green’s  function  based 
upon  the  Helmholtz  operator.  Solutions  of  this  in¬ 
tegral  equation  proved  to  be,  on  the  average,  of  the 
same  order  of  magnitude  as  the  CFD  -  integral  methods 
but  resulted  in  a  smoother  error  curve  over  the  calcu¬ 
lated  values.  Also  the  Green's  integral  calculations 
proved  quicker  than  either  the  integral  equation  method 
or  using  the  CFD  quotient. 


Case  III 

The  third  and  final  problem  investigated  in  one- 
dimension  was  the  inhomogeneous  Cauchy-equation  with 
homogeneous  Dirichlet  boundary  conditions.  Expressed 
mathematically,  the  problem  is 


(208) 


with  boundary  conditions 


(209) 


Both  previous  results,  that  is,  converting  the 
differential  equation  into  integral  form  and  approx¬ 
imating  the  solution  through  numerical  integration 
and  using  the  method  of  finite  differences  have  proven 
to  be  equivalent.  One  possible  reason  for  this  could 


» 

1 
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be  due  to  the  linear  kernel  which  has  resulted  in  both 
integral  conversions.  Equation  (208)  was  found  not 
to  have  a  linear  kernel  upon  integral  conversion. 

The  properties  of  this  equation  are  now  investigated. 


Analytical  Solution 

The  complete  general  solution  to  (208)  is  found 

to  be 

/)  .  Bxa  +  A?  (210) 

3T  10 


Using  boundary  conditions  (209) »  the  constants  /J  and 
B  are  evaluated.  The  unique  solution  is 


(211) 


and  is  graphed  in  Figure  14.  Note  that  the  solution 
is  asymmetric. 


Integral  Solution 

The  Cauchy  differential  equation  can  be  expressed 
in  integral  form  by  first  rearranging  terms  and  solving 
for  the  second  derivative.  Equation  (208)  becomes 

(212) 

Now,  integrate  both  sides  with  respect  to  x  over  the 
interval  (l,x).  Equation  (212)  becomes 

J(x'-e*£>)dx  (213) 
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Integrating 


where  the  constant  C  represents  . 

d* 

a  second  time  over  the  same  limits  makes  (213)  equi¬ 


valent  to 


y*) 


*  C  J"<U  + 


(214) 


By  using  the  results  of  (6)  for  the  multiple  integral 
and  performing  the  integration  of  the  first  integral, 


equation  (214-)  becomes 


rx 

y(K);Cx-C  +  Jrx-O]?1  ^ 


(215) 


The  constant  C  can  be  evaluated  by  using  the  boundary 


condition 


tjte.)  -  o  -  2c  v 


(216) 


Solving  for  C  yields: 


c  =  i 


Jcs-sj/'s*  +  £i<aj  Jj 


(217) 


Using  this  expression  for  C,  the  complete  solution 
for  yl*).  given  by  equation  (215)  becomes 

°U  (218) 

i  1 

Equation  (218)  can  be  put  in  the  simplified  form 

r3 

(219) 


by  using  the  abbreviation 


83 


AD-AOSO  *18  UR  FORCE  INST  OF  TECH  NRISHT-FATTERSON  A  FI  OH  SCHOO — ETC  F/a  It/1 
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K<x,5)= 


a 

(X-3Xf->) 

a 


(220) 


Equation  (  219)  represents  the  equivalent  integral  form 
of  the  one-dimensional  Cauchy  equation.  The  integral 
equation  can  be  rearranged  by  multiplying  through  by 
and  seperating  tne  integral.  Combining  terms,  equa¬ 


tion  (  219)  is  also  equivalent  to  the  expression 

3  3 

y(x;  -  a 


(221) 


The  integral  on  the  right  simplifies  by  replacing  the 
kernel  with  its  values  over  the  appropriate  limits. 


The  integral  becomes 


3 

f(y-)Q  -3)  {*■ 


(222) 


Integrating  both  integrals  over  §  ,  equation  ( 222) 


becomes 


J<_  _  U>  X  v  13 

»2  3  ^ 

Equation  ( 221)  then  becomes 

**  =  £  -  '#*  +  ^ 


(223) 


(224) 


It  is  now  possible  to  use  the  trapezoid  approximation 
for  the  left  side  integral  and  equation  (  224)  becomes 
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-  £  £  fl  ***»■.)  = 


n  *  'fx*f  (.225) 


where  the  points  y^  have  been  chosen  over  the  in¬ 
terval  (1,3).  In  matrix  notation  expression  (225) 
transforms  into  the  system  of  linear  equations 


fl-  -  x 


(226) 


where  D  corresponds  to  the  weighing  coefficients 
used  in  the  numerical  integration,  ft  *  is  the  modified 
kernel  with  denominator,  and  X  corresponds 

to  the  function  vector  defined  by  (  223) •  Equation  (  226) 
can  now  be  solved  using  Gaussian  elimination. 


Green's  Integral 

In  addition  to  the  integral  representation  it  is 
also  possible  to  express  the  solution  of  •j1*)  by  the 
Green’s  integral  equation.  From  Appendix  C,  the  Green's 
function  for  the  one-dimensional  Cauchy  operator  was 


found  to  be 


-*■) 

Tt(s  ~?)(x  ~  &)  x>f 


(227) 


Therefore  the  equivalent  expression  for 
expressed  as 


y*> 


can  be 
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(228) 


where  (JooO  is  as  defined  above.  Using  the  trapezoid 
rule,  equation  (228)  becomes  equivalent  to  the  quantity 


6/  -  £  dk  * 

K-  l 


t 

K 


(229) 


or  in  matrix  notation  |j(*)  appears  as 


(230) 


where  D  and  0  are  both  matrix  vectors  of  dimension 

— •  a 

N  X  N  and  X  is  the  column  vector  of  dimension  N  X  1 . 


Finite  Differences 

The  method  of  finite  differences  was  used  to  com¬ 
pare  numerical  differentiation  with  the  integral  and 
Green’s  methods  outlined  above.  Again,  the  derived 
central  finite  difference  quotient,  equation  (97),  was 
substituted  into  the  one-dimensional  Cauchy  equation. 
After  collecting  terms,  equation  (208)  became  equiva¬ 
lent  to  the  expression 

+  j/;.,  =  xfh*  (23D 

The  matrix  has  the  tridiagonal  form 


( 
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Using  the  boundary  conditions,  O  and  0  *  ex¬ 

pression  (232)  can  be  rewritten  in  the  form 
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The  dimensions  of  all  three  matrices  have  been  re¬ 
duced  to  N  -  2.  In  concise  matrix  form  expression 
(233)  can  be  written 

|cj  -  R  (234) 

where  &  corresponds  to  the  given  band  coefficient 
matrix  and  R  represents  the  column  function  vector 
defined  by  the  right-hand  quantity  in  (233  )  • 

Relative  Error 

For  a  comparison  of  the  relative  errors  asso¬ 
ciated  with  the  three  numerical  methods  see  Figure  15* 

The  relative  error  percentages  are  tabulated  in  Table  VIII . 

Table  VIII.  Relative  Error  {%)  -  Computer  Solution 
for  x*y"-2y=x$ ,x=  2 


Number  of  Iterations 
per  Interval 

Integral 

Equation 

Green's  Integral 
Method 

CFD 

Method 

2 

2.402402 

.150144 

6.306306 

4 

.732556 

.009369 

1. 718748 

6 

.3^4530 

.001815 

.784274 

8 

.198287 

.000586 

.445971 

50 

.005241 

.000000 

.011593 

Error  Analysis 

Because  both  the 

integral 

and  Green’s  represents- 

tions  are  equivalent  either  expression  can  be  used  to 
predict  an  error  bound  using  the  trapezoid  rule.  The 
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Comparison  of  Relative  Errors 


Green's  integral  representation  is  used  once  again 
because  it  does  not  involve  the  function  yff)  ,  making 
the  error  calculations  easier  to  work  with.  The  Green’s 
integral,  equation  (228),  can  be  written 

I 

+  J«j  <235) 

x 

where  G(*,e)  has  been  substituted  over  the  appropriate 
limits.  The  error  will  be  computed  at  the  point  X**2. 
and  expression  (235)  simplifies  to 

y<*)  •  -  ^ -U-^JV-Ods  -  v.-wws  jf f -i)  djj-  (236) 

The  maximum  error  for  both  integrals  is  the  sum  of 
errors  associated  from  each.  To  find  the  upper  bound 


on  the  errors  all  errors  are  assumed  to  be  additive. 
The  upper  bound  can  now  be  associated  with  fit). 
where  ■fio  corresponds  to  the  sum  of  the  two  inte¬ 


grands  above: 


/<$)*  +  >■*"*'(£-§)  (237> 

To  calculate  ^  the  second  derivative  of  fti)  is 
needed.  This  is  given  by 

*.$'  (238) 

The  value  of  fit)  at  the  point  fsZ  is  found  to  be 
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1 


If}?)  I  =  H 


(239) 


and  the  associated  error  is  of  the  order 

¥  (k-a-)A^Z  _  .(.UH  (240) 

\Z. 

This  error  represents  only  a  crude  approximation  to 
those  errors  associated  using  trapezoid  integration 
and  can  be  misleading.  All  errors  in  arriving  at 
(240)  were  assumed  to  be  additive.  Because  many  of 
these  errors  cancel  the  predicted  values  will  always 
be  larger  than  observed,  as  seen  in  Table  IX. 

The  relative  error  was  also  examined  at  each  array 
element  when  the  interval  (1,3)  was  divided  into  50 
partitions.  The  errors  associated  with  the  integral 
solution  are  graphed  in  Figure  16.  The  Green’s 
integral  solution  converged  to  the  exact  solution  at 
every  array  point,  to  six  decimal  place  accuracy  and 
was  not  graphed.  The  relative  error  associated  with 
the  CFD  method  was  also  examined  at  each  array  element 
and  is  graphed  in  Figure  17.  Note  that  both  sets  of 
data  display  similar  trends;  the  higher  numbered  array 
elements  have  lower  errors  as  compared  to  the  lower 
numbered  array  elements.  This  trend  can  be  explained 
by  examining  the  method  of  solution.  Both  methods 
require  the  simultaneous  solution  of  a  set  of  N  alge¬ 
braic  equations,  49  equations  for  the  CFD  method  and 
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Relative  Errors  of  Array  elements 
Integral  Equation  Method 


51  for  the  integral  approach.  Both  sets  are  solved  by 
Gaussian  elimination.  As  was  covered  earlier  the  last 
array  elements  are  evaluated  first  and  therefore  have 
the  smallest  round-off  errors.  As  the  elements  above 
are  computed  from  the  know  elements  below  there  is  an 
accumulation  of  air  with  each  new  row  evaluated.  This 
accounts  for  the  trend  in  increasing  error  as  the  last 
elements  are  calculated,  corresponding  to  the  lower 
numbered  array  elements. 


Simpson's  Comparison 

The  trapezoid  rule  was  compared  to  the  Simpson's 
rule  for  the  same  array  elements  as  above,  see  Figures 
18  and  19.  Note  the  downward  trend  for  the  relative 
error  values  at  the  higher  indexed  array  elements  in 
both  figures.  Also,  the  same  sinusoidal  character  is 
exhibited  over  the  odd  -  even  partitions  as  was  observed 
in  Figure  8.  In  addition,  Simpson's  rule  gives  a  better 
approximation  to  the  calculated  values  at  the  odd 
numbered  intervals  than  the  trapezoid  rule.  This 
results  because  the  kernel  associated  with  the  integral 
equation  is  no  longer  a  linear  function,  but  rather 
can  be  expressed  as  the  modified  kernel 


fx-Oft-.?) 

ti3- 


X  <  § 

*>S 


(241) 
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50PTS.  TRAPEZOID  RLE 


SOFTS.  SIMPSON  RULE 


Quadratic  interpolation,  i.e.,  parabolas,  provide  a 
more  accurate  estimate  of  the  non-linear  kernel  than 
using  trapezoids. 

Times  of  Solution 

The  two  numerical  integration  routines  were  com¬ 
pared  to  the  GFD  method  for  speed.  The  results  of 
this  comparison  are  given  in  Table  X. 

Table  X.  Computer  Time  -  x  y"-2y=x^ 


Method 

Compilation  Time 
( sec) 

Execution  Time 
( sec) 

Usage 

(Kilo-word  sec) 

Integral 

Equation 

0.799 

1.186 

194.087 

Green's 

Integral 

0.662 

0.847 

117.473 

CFD 

0.566 

0.557 

169.940 

The  CFD  method  was  executed  2.1  times  faster  than 
the  integral  equation  method  and  1.5  times  as  fast  as 
the  Green’s  integral  approach.  In  addition,  the  CFD 
method  was  compiled  29.2$  faster  and  used  12.4fS  less 
central  storage  than  the  integral  method.  Also,  the 
CFD  method  was  complied  14.5$  faster  and  used  30*8$ 
more  central  storage  compared  to  the  Green’s  integral 
equation. 

Conclusions 

The  one  significant  feature  about  the  one-dimen- 
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sional  Cauchy-type  equation  is  the  fact  that  its  so¬ 
lution  is  asymmetric.  Converting  the  original  differ¬ 
ential  equation  into  integral  form  still  results  in  a 
linear  weighted  Green's  kernel  but  multiplied  by  a 
non-linear  function  of  $  .  Also  the  Green's  integral 
equation,  for  the  differential  operator,  has  a  non-sy- 
mmetric  Green's  kernel.  Results  by  the  Green’s  integral 
equation  proved  superior  to  the  other  two  methods. 

The  integral  equation  method  also  gave  more  accurate 
solutions  as  compared  to  the  CFD  method. 


IV.  Two  Dimensional  Numerical  Methods 


The  Steady-State  Heat  Conduction  Problem 

The  last  case  investigated  in  this  thesis  will  be 
the  two-dimensional  steady- state  heat  conduction  prob¬ 
lem  with  inhomogeneous  Dirichlet  boundary  conditions. 
The  internal  temperature  within  the  rectangular  plate, 
shown  in  Figure  20,  staisfies  the  partial  differential 


equation 


Vr*T(*,u)  s  -FCx.. 


(242) 


where  TCv.y)  is  the  temperature  at  the  point 
and  represents  the  internal  heat  generation. 


Fig.  20  Two-Dimensional  Steady-State 
Heat  Conduction  Problem 


The  boundary  conditions,  from  Figure  20  are 

T(\o)  -  T(0,jj)  S  T(x,b)  =  Ot  T(a,V|)iTc 


(243) 
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where  To  represents  some  constant  temperature  along 
the  edge  OL  .  For  all  problems  considered,  F6c,j)  is 


assumed  constant  over  the  plate  and  has  the  value 


Rx,^)  =  i 


(244) 


Analytical  Solution 

The  analytical  closed-form  solution  of  equation 
(242),  Poisson's  equation,  can  be  expressed  as  the  sum 
of  the  two  harmonic  functions;  the  first  satisfies 
equation  (242)  with  homogeneous  Dirichlet  boundary 
conditions,  and  the  second  function  satisfies  the 
Laplacian  of  (242)  with  the  inhomogeneous  boundary 
condition  at  T(a,3)  .  The  complete  solution,  using 
condition  (243)  and  satisfying  the  given  boundary 
conditions,  is  derived  in  Appendix  D  and  found  to  be 

+  5  i  ±  H  V 5,n  t#'4  T)  (245) 

Note  that  the  solution  is  a  double  Fourier  series, 
a  summation  process  which,  if  slowly  convergent  will 
take  considerable  time.  For  the  problem  under  inves¬ 
tigation  the  temperature  at  16  interior  nodal  points 
waS  examined.  Also,  instead  of  using  the  rectangular 
area  the  analytical  model  was  altered  into  a  square 
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of  unit  dimension.  From  equation  (245)  the  internal 
temperature  was  computed  at  81  interior  grid  points, 
a  computation  requiring  70.9  seconds  for  six  place 
decimal  accuracy. 

Finite  Differences 

The  solution  to  the  two-dimensional  steady-state 
heat  conduction  problem  by  finite  difference  methods 
is  simply  an  extension  of  the  technique  used  in  the 
one-dimensional  cases.  The  first  step  is  to  set  up  a 
grid  system  of  interior  points  within  the  plate,  shown 
in  Figure  21 . 


Fig.  21  Interior  Grid  Network 


The  Laplacian  of  equation  (242)  can  be  approximated 
by  expressing  the  two  second  derivatives  of  the  temp¬ 
erature  as  a  difference  equation.  The  second  deriva¬ 
tives,  using  equation  (97)  and  the  interior  nodal 
system  in  Figure  21  become 


f 


(246) 


£«  **b>  =  T(x-K^)  -  £T(x,»j) 

jx*  ^  - 


and 

£t*x,ij)  TCx^-fe) -2T(x,^)  +T(xlu+h) 

*r  - if - 


(247) 


Substituting  these  approximations  back  into  the  heat 
equation, ( 242) ,  and  the  combining  terms  gives  the 
expression 


T(x-V),<j)  +  T(x+V>,<j) 


+  T(*.^-W)  +  T(x,^V,)  =  (2^8) 


I 

i 


Once  again  the  finite  difference . method  results  in  a 
set  of  linear  algebraic  equations. 
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Fig.  22  Interior  Nodal  Arrangement 
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Using  the  nodal  point  coordinate  system,  shown  in 
Figure  22,  and  the  boundary  conditions  given  by  ex¬ 
pression  (243),  the  following  matrix  representation 
of  equation  (248)  is  obtained 
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where  T*,*,  corresponds  to  the  nodal  temperatures 


shown  in  Figure  22,  and  corresponds  to  the  heat 
generation  function  at  the  nodal  points  (n,m). 


81  Interior  Nodes 

In  addition  to  the  4x4  grid  system  for  the  16 
interior  nodal  points,  the  finite  difference  method 
was  used  to  examine  the  results  when  the  spacing 
between  adjacent  nodes  decreased.  The  original  4x4 
grid  was  subdivided  in  half,  and  the  temperature  at 
81  nodes  was  in  a  new  9x9  grid  network.  Instead  of 
generating  81  algebraic  equations  note  that  there  is 
a  pattern  in  the  diagonal  vectors  appearing  in  the 
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coefficient  matrix  of  expression  (249).  If  the  sequen¬ 
tial  location  of  the  zero  elements  along  any  diagonal 
is  known  then  it  is  a  relative  easy  manner  to  generate 
the  entire  diagonal  vector.  Also  because  there  are  at 
most  five  row  elements,  from  equation  (248),  in  any 
row  of  the  coefficient  matrix  only  five  diagonal  vec¬ 
tors  have  to  be  stored.  A  comparison  of  the  two  re¬ 
sults  follow. 

Results 

The  solution  of  (249)  by  Gaussian  elimination, 
using  the  method  of  finite  differences  is  compared  to 
the  exact  analytical  solution  in  Table  XI .  Due  to  the 
symmetry  in  the  original  problem  only  half  the  nodes 
are  tabulated. 

Table  XI.  Finite  Difference  Solution  for 
Square,  with  Heat  Generation 


Node 

Exact 

Solution 

Finite  Diff 
16  Interior 

Error 

Finite  Diff  Error 

81  Interior 

2,2 

4.331240 

4.578788 

-.247548 

4.446669  - 

.115429 

2,3 

6.969133 

7.243636 

-.247503 

7.115880  - 

.14674 7 

3,2 

10.555374 

11.031515 

-.476141 

10.755899  - 

.200525 

3,3 

16.772428 

17.112121 

-.339693 

16.969429  - 

.197001 

^,2 

21.727495 

22.395152 

-.667657 

22.012913  - 

.285418 

4,3 

33.090082 

33.021212 

-.068870 

33.166982  - 

.076900 

5,2 

45.599465 

45.487879 

+.111586 

45.621920  - 

.022455 

5,3 

60.555373 

59.516364  +1.039009 

60.306397  + 

.248976 

Comment  -  Finite  Differences 


The  method  of  finite  differences  in  two-dimen¬ 
sions  can  be  viewed  as  an  averaging  process,  using 
the  known  boundary  conditions  to  generate  values 
at  other  interior  points.  The  best  results  are  ob¬ 
tained  when  all  nodal  points  are  symmetrical  then 
the  boundary  values  at  the  edges  are  'weighed'  more 
evenly  at  the  interior  points.  To  demonstrate  this 
consider  the  heat  conduction  problem  as  before;  the 
internal  temperatures  at  the  point  (.5* *5)  will  be 
computed  by  two  seperate  schemes.  The  first  pro¬ 
cedure,  shown  in  Figure  23 »  is  to  use  three  adjacent 
points  located  parallel  to  the  x-axis  and  spaced 
\l=  apart.  The  temperature,  T,  ,  calculated  by 
the  difference  equation  (248)  is 

-J-  s  23.5  <»SS?2 


Fig.  23  Finite  Difference  Scheme 

Using  Three  Points  Parallel 
to  the  /-Axis 


(250) 
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The  procedure  shown  in  Figure  24,  is  to  use  three  ad¬ 
jacent  points  again,  but  located  parallel  to  the  y-axis, 
spaced  apart.  T,  calculated  by  the  difference 

equation  (248)  is  now  found  to  be 

T;  *  at.sjm 5  (25D 


From  the  analytical  solution,  summed  over  10,000  sep¬ 
arate  iterations,  the  value  is  24.926328.  The  temper¬ 
ature  is  observed  to  be  greater  where  the  boundary 
condition  7o  is  used  three  times.  It  is  interesting 
to  note  bhat  using  only  one  interior  node  and  equation 
(  248) ,  T,  is 


T  =  ZS.OLZS  (252) 

'I 

This  value  is  only  .55#  in  error  from  analytical  value 
at  (.5, .5). 


Fig.  24  Finite  Difference  Scheme 

Using  Three  Points  Parallel 
to  the  y-Axis 
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Introduction  -  Splines 


The  last  method  to  be  presented  in  two-dimen¬ 
sions  is  the  integral  solution  to  the  steady-state 
heat  conduction  problem  using  the  method  of  cubic 
splines. 

The  integral  solution  for  heat  conduction  over 
a  rectangular  region,  with  homogeneous  Dirichlet 

boundary  conditions,  can  be  expressed  in  the  form 

a  b 


^  Jvj  (253) 

where  G(m  ;^,n)  is  the  two-dimensional  Green’s 
function  for  Laplacian  (Ref  13s  520-523).  If  the 
exact  Green's  function  is  known,  and  is  sep- 

erable,  equation  (253)  can  be  evaluated  using  the 
one-dimensional  techniques  already  considered.  Many 
times  though,  for  problems  other  than  rectangular 
symmetry,  the  Green's  function  is  extremely  diffi¬ 
cult  to  compute;  the  Green’s  function  depends  upon 
the  boundary  symmetry,  for  complicated  regions  Gfai'yn) 
can  be  exceedingly  complex.  An  alternate  integral 
form  of  solution  can  be  expressed  for  T<*'3)  using  splines. 
The  method  to  be  presented  was  first  proposed  by 
Hadjin  and  Krajcinovic  using  cubic  splines  in  integral 
form  for  solving  elliptical  partial  differential 
equations  (Ref  2«  513-539). 
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Spline  Functions 


Because  the  proposed  method  uses  cubic  splines 
a  brief  introduction  into  spline  theory  is  necessary. 
Spline  functions  are  simply  polynominals  used  to 
approximate  a  given  function,  ■f(x)  ,  over  an  interval 
(a,b).  However,  instead  of  using  a  single  polynominal 
over  the  entire  domain  the  interval  (a,b)  is  divided 
into  subintervals, 


*,  <■  x2  •  •  •  <  xN  =  \>  (25*0 

with  a  different  polynominal  representing  •f(x)  over 
each  subdivision.  By  definition,  a  cubic  spline  is  a 
continuous  piecewise  function  having  continuous  first 
and  second  derivatives  everywhere  on  the  interval  (a,b) 
and  is  represented  by  a  polynominal  of  degree  three  or 
less.  Hence,  the  spline  Vo  consists  of  cubic  poly¬ 
nominals,  one  in  each  of  the  subintervals  (*k-,,Xk)  . 

The  cubic  spline  V*>  ,  which  approximates  "foe)  in 
each  interval  (**»**)  can  be  uniquely  determined  by  in¬ 
sisting  that  the  value  of  the  spline  and  its  first 
derivatives  are  equal  to  the  value  of  foO  and  its  first 
derivatives  at  each  node,  •. 


Sta)  =  f  fXi)  *n«j  sfao  -  f(Xi)  ( is  *‘l>  K) 


(255) 


In  addition,  to  insure  a  smooth  fit  of  the  spline 
function  across  the  interval  (a,b),  the  values  of  the 


109 


first  and  second  derivatives,  on  either  side  of  the 
nodal  points  will  be  assumed  equal,  that  is 

and 


(256) 


S^(v)  r  S  /xkO 


(257) 


Using  the  conditions  required  by  expressions  ( 255) » 

(  256),  and  (  257),  the  cubic  spline  representing  S Oc) 
can  be  shown,  Appendix  E,  to  be  equivalent  to 

S..(x)  =  /'V.  (<*-*)’  + 

K  <-*>«  ~ ChT 

,  xjbK-a.nl  - ^(>,3.,,-  J.) -fc.M,,, - „y^/4j 

where  the  symbols  are 

h*  *  Xfc-X*-, 

ito*  f  uh) 

(Ref  14}  296-298) 

By  using  the  continuity  requirements  for  the 
spline,  from  equation  (  2 55),  the  values  of  can  be 
determined,  namely  it  is  required  that 


(258) 


V£M*-»  t  Ivt  VfeaMfc  t  ViKt,  Mm 


nr 


r 


Expression  (259)  represents  a  system  of  N-l  linear 
algebraic  equations  with  unknowns  M*»  A^,  By 


(259) 


110 


assigning  the  values  of  zero  to  ^and  V  because 
the  physical  sline  can  be  assumed  to  be  straight 
outside  the  interval  (a,b),  i .  e .  ,t(X)=  O  for  X<^ 
or  X>V  ,  the  other  values  can  be  de¬ 

termined.  Using  the  spline  representation,  equation 
(  258),  and  the  continuity  requirements,  expression 
(  259) ,  the  method  of  splines  in  the  solution  of  integral 
equations  can  now  be  introduced. 

Integral  Method  of  Solution 

The  steady-state  heat  conduction  problem  has 
the  form 

iTr*,)  +  i'T/K.'j)  «  -  (260) 

In  the  method  proposed  by  Hajdin  and  Krajcinovic  the 
highest  derivatives  of  the  function  are  chosen  as 
the  unknown,  that  is  let 

oni  (261) 

1) 

Substituting  equations  (261)  back  into  equation  (260) 
yields 

~  (262) 

Differential  equations  can  now  be  written  equations 
(261)  for  a  fixed  value  of  X*  X*  and/or  *j«  , 

£%*'$*)  * ^  (263) 

«X  in1 


ill 
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Equations  (263  )  can  now  be  converted  into  integral 
form  using  the  Green's  function  for  the  differential 
operators  in  equations  (263  ) .  The  equivalent  integral 
representation  of  equations  (263  )  are 

CL 


T(x.>)=  +  Ikl 


(264) 


and 

In  (265) 

o 

where  G«,S)  and  are  the  one-dimensional  Green's 

functions  for  the  second  order  differential  operators 
appearing  in  equations  (263  ) .  Note  that  the  partial 
differential  problem  has  been  converted,  along  a  con¬ 
stant  line  K-  X*  and/or  y*  JJn  »  into  two  one-dimen¬ 
sional  integral  equations.  The  unknowns  pis,*.)  and 
V’lfc’i)  are  now  approximated  using  cubic  splines,  given 
by  equation  (258  )  and  integrated  in  equations  (264  ) 
and  (265  ) . 


Spline  Integration 

Using  the  spline  representation,  equation  (264) 
c  ome s 

p  ♦(q-KQy  O.T  t  T*xk  (266) 
*  l  H5  £4  3  *  J  ~  J 

and  equation  (  265)  becomes 
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-fc-a-fl<n«x -(«-%)>.  i>-»1r„4  (267 

(  MS  2M  3  Z  j 


where  =  length  of  a  side  of  the  square  plate 

=  location  within  region 

=  value  of  the  second  derivative  of  the 
spline  at  (*K|yn)  . 

\‘%= 

Two  expressions  for  the  temperature  now  exist  at  the 
point  (Xk.tjq)  given  by  expressions  (  266)  and  (  267). 
Because  the  temperatures  must  be  equivalent  at  these 
locations,  equations  (266)  and  (267)  can  be  equated. 
Using  the  continuity  requirements  for  the  second 
derivatives  of  the  spline  functions,  given  by  ex¬ 
pression  ( 259) ,  both  /MK  and  ^»|  can  be  replaced  by  a 
multiple  of  p^  and  .  In  addition,  a  second  condi¬ 
tion  exists  between  and  ^  given  by  equation  ( 262) . 
Therefore,  using  the  two  algebraic  equations,  (266) 
and  (267),  and  equation  (  259) ,  the  unknowns  and 
can  be  evaluated.  Substituting  these  values  for  p^ 
or  back  into  equations  (266)  or  (267),  along  with 
the  substituted  values  of  M*  or  will  yield  the 
temperature  at  the  point 


Results 

The  computed  values  at  the  same  interior  nodes 
used  in  the  method  of  finite  differences  is  given  in 
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Table  XII.  Observe  how  the  error  values  increase 
as  the  node  points  move  away  from  the  center 

Table  XII.  Cubic  Spline  Solution  for 

Square,  with  Heat  Generation 


Node 

Exact 

Solution 

Cubic  Spline  Error 

16  Interior 

2,2 

4.331240 

10.040000 

-5.708760 

2,4 

6.969133 

10.617943 

-3.648810 

3,2 

10.555374 

18.903657 

-8.348283 

3,3 

16.772428 

20.049333 

-3.276905 

4,2 

21.72749 5 

28.332229 

-6.604734 

4,3 

33.090082 

30.049333 

+3 . 040749 

5,2 

45.599465 

40.044000 

+5.555465 

5,3 

60.555373 

42.332229 

+18.223144 

Besides  the  16  interior  points,  the  grid  spacing  was 
halved  to  observe  the  effects.  Instead  of  changing 
the  computed  values  at  the  node  points  the  exact 
same  values  as  before  were  calculated,  demonstrating 
that  unlike  the  finite  difference  method,  the  method 
of  cubic  splines  does  not  improve  by  reducing  the 
grid  spacing.  It  is  interesting  to  note  that  for  the 
most  symmetrically  located  point,  (.5*  *5)*  "the  cal¬ 
culated  temperature  is  within  .50%  of  the  actual 
analytical  solution. 

Conclusions 

The  method  of  finite  differences  in  two-dimen¬ 
sions  gives  excellent  results,  the  largest  relative 


error  was  5*72#  at  the  node  2,2.  Halving  the  interval 
spacing  resulted  in  a  decrease  of  this  error  to 
2.67 

The  integral  method  using  cubic  splines  appears 
to  work  well,  but  only  at  the  most  symmetrically  lo¬ 
cated  point.  The  major  fault  with  this  method  comes 
when  converting  the  two-dimensional  problem  into  one¬ 
dimensional  form.  Much  of  the  physical  significance, 
the  edge  effects  from  the  boundaries  and  the  inhomo¬ 
geneous  boundary  conditions,  have  been  almost  com¬ 
pletely  neglected.  For  the  particular  part  of  the 
solution,  given  by  equation  (266),  no  account  for  the 
value  is  used  though  in  reality  the  particular 

solution  is  both  a  function  of  two  variables,  X  and  vj  . 
Also  there  does  not  appear  sufficient  coupling  to 
assume  that  the  two  seperate  splines  and  condition 
(262)  should  result  in  the  same  temperature  at  the 
point  »  equation  (267)  is  based  upon  homogeneous 

boundary  conditions  and  equation  (266),  for  the  same 
point,  incorporates  two  distinctly  different  boundary 
conditions.  In  contrast  to  what  the  authors  Hajdin 
and  Krajcinovic  claim,  i.e.,  that  problems  of  potential 
theory  may  be  easily  and  conveniently  solved  employ¬ 
ing  the  method  of  cubic  splines  and  integral  con¬ 
version,  the  results  that  were  found  do  not  justify, 
for  the  reasons  cited,  the  authors  claims. 
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V.  Conclusions  and  Recommendations 


Conclusions  One-Dimensional  Cases 

The  purpose  of  this  thesis  has  been  to  inves¬ 
tigate  and  compare  numerical  integration  techniques 
with  methods  based  upon  numerical  differentiation, 
the  finite  difference  method.  In  one-dimension  the 
results  of  this  study  have  been  somewhat  inconclusive; 
in  the  three  examples  investigated  the  first  two 
have  proven  to  yield  equivalent  solutions  for  both 
the  integral  equation  approach  and  the  method  of 
central  finite  differences.  This  is  due  to  the  simi¬ 
larity  of  the  two  examples,  the  Helmholtz  equation  in 
Case  II  was  but  a  modified  version  of  the  particular 
equation  for  Case  I.  In  addition,  the  Green's  inte¬ 
gral  approach  for  Case  II  was  found  to  yield  solutions 
on  the  same  order  of  error  as  the  integral  and  CFD 
methods.  In  contrast  however,  using  the  particular 
Cauchy- type  equation  for  Case  III,  three  entirely  diff¬ 
erent  results  for  each  method  was  observed.  The  Green's 
integral  approach  proved  to  be  the  most  accurate,  for 
only  one  iteration  over  the  interval  (1,3)  the  com¬ 
puted  value  was  found  in  error  by  only  .156$.  The 
integral  method  however  for  the  same  equation  gen¬ 
erated  approximately  15  times  the  error  of  the  Green’s 
integral  but  for  a  single  iteration  the  error  was 
still  less  than  3%.  The  CFD  method  had  more  than 
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twice  the  error  of  the  integral  method,  6.3 %• 


Another  similarity  between  Case  I  and  Case  II 
was  the  fact  that  bcth,  upon  integral  conversion,  had 

linear  kernels.  This  can  be  attributed  to  the  second 

J* 

order  differential  operator,  ,  appearing  in 

both  differential  equations.  Case  III  however,  had 
a  non- symmetrical  kernel  and  asymmetric  solution,  in 
contrast  to  the  first  two  cases. 

Conclusions  Two-Dimensional  Cases 

The  method  of  finite  differences  in  two-dimen¬ 
sions  gives  excellent  results,  at  all  points,  compared 
to  the  analytical  solution  of  the  steady-state  heat 
conduction  problem.  In  contrast,  the  integral  approach 
using  cubic  splines  falls  short  of  its  expectations  at 
all  interior  points  except  the  center.  The  major 
drawback  of  the  spline  technique  occurs  when  the 
original  two-dimensional  problem  is  expressed  as  two 
separate  one-dimensional  cases.  The  coupling  scheme 
between  the  splines  is  no  longer  valid  due  to  diff¬ 
erent  Dirichlet  boundary  conditions  existing  at  the 
plate  edges.  Another  serious  drawback  to  the  spline 
method  is  that  it  does  not  improve  as  the  grid  net¬ 
work  is  decreased.  In  contrast,  the  finite  difference 
solutions  improve  as  the  grid  spacing  is  reduced. 


( 
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Recommendations 


More  study  needs  be  given  to  a  wider  variety  of 
one-dimensional  problems,  using  both  linear  and  non¬ 
linear  kernels.  Higher  order  approximations  to  the 
derived  difference  quotients  should  also  be  included 
to  reduce  truncation  errors  in  the  finite  difference 
technique.  In  addition,  other  integration  schemes 
could  also  be  used  to  see  how  well  specific  integral 
techniques  compare  to  the  numerical  differentiation 
techniques.  For  the  two-dimensional  analysis  a  two- 
dimensional  spline  function  could  be  tried.  The 
disadvantage  of  working  in  two-dimensions  though  is 
that  the  exact  form  of  the  Green's  function  must  be 
known,  for  complicated  geometries  this  could  prove 
extremely  difficult. 

In  conclusion,  there  are  many  areas  in  which 
further  studies  may  lead  to  more  fruitful  results. 
This  thesis  has  but  presented  a  brief  analysis  of 
several  numerical  methods  for  the  solution  of  diff¬ 
erential  equations  in  both  one  and  two-dimensions. 
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Appendix  A 


Kernel  Properties  of  the  Equation  y  +  ajg*  b a-o 
for  al 

A  differential  equation  of  the  form 

t  A4*  *  8u  = 
dx2  ^ 

with  homogeneous  Dirichlet  boundary  conditions 


(268) 


y<o)  = 


y(t)  - 


(269) 


will  have  a  non- symmetrical  kernel  when  expressed  in 
integral  form  if  4*o  .  In  order  to  demonstrate 
this,  rewrite  equation  (268)  in  the  form 


fe  =  "Afe  -B: 


(270) 


and  integrate  both  sides  over  the  interval 
Equation  (270)  now  becomes 


h'x 

<)x 


(271) 


or 


l 

jijjx)  -  4y(o)  _  -tfyoo  (272) 

o 

Let  C  donote  the  constant  4f  and  integrate  equation 
(272)  a  second  time  over  the  same  limits  as  before. 
Expression  (272)  becomes 


Jf'  X  X  *  * 

n  C>  0  00 


(273) 


121 


Now,  using  the  results  of  equation  (6  )  and  changing 
the  dummy  variable  of  integration  in  the  second  in¬ 
tegral  on  the  right-hand  side  of  equation  (273  )  it  is 
found  that 

X  K 

vj(x)  =  Cx  - 

o  o 

Applying  the  boundary  conditions  at  evaluates 

the  constant  C,  that  is 

|  1 

y(j)=r  O  «  -  Aj^o4§  ~ 

D  o 

c=  ty '»*•«-«>} 

and  equation  (274),  using  C  above  becomes 

i  x 

Expressing  the  first  integral  as  |  r  f  t  [  becomes 

x  j.  *  •  * 

o  J  X 

The  kernel  in  the  above  integral  equation  expression 
can  be  identified  by  making  the  substitution 

where  the  abbreviation  for  the  kernel  is 


(274) 


(275) 


(276) 


(277) 


(278) 


(279) 


K  {*,$)  =. 


J-  ■»  Bf  u-f) 

A-f±  -  »i(f-0 


(280) 


r 


Note  that  the  kernel  is  both  non- symmetric  and  dis¬ 
continuous  at  X*5,  unless  -  O 
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Appendix  B 


Determining  the  Green's  Function  for  the 
One-Dimensional  Helmholtz  Equation 

The  differential  equation  for  the  one-dimensional 


Helmholtz  equation  is 

&  +  >rr  -f(x) 


(281) 


The  Green's  function  for  the  differential  operator, 
i.e.,  the  left  side  of  expression  (281),  must  sat¬ 
isfy  the  following  differential  equation 


-  -  $(*-%) 


(282) 


Because  G6c,$)  can  be  thought  of  as  a  function  of  only 
one  variable  for  a  constant  f  ,  assume  that  the  Green’s 
function  can  be  expressed  by  Gt(x)  and  Ga<*>  over  some 
interval  (»,*)  .  The  Green’s  function  will  be  of  the  form 

GU,«JG'(X)  K<*  (283) 

-  G2(X)  x>s 

for  G|(x)and  Ga(x)  will  each  satisfy  the  homogeneous 

equation  of  expression  (282  )»  that  is 


dXl 


(284) 


^%(x)  +  X*  ^Cx)  ST  O 


(285) 


The  solution  of  the  homogeneous  equations  (284  )  and  (285  ) 
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yield,  for  G,  (xl  and  Ga(x)  ,  the  Green's  function 


G(x,s)  - 


A  cos  Xx  -V  Bsm  x  x 
DcesXx  t  Esm^x 


*<5 


(286) 


The  Green's  function  has  the  property  of  being  zero 
at  the  end  points  and  therefore  two  of  the  constants 
in  expression  (286)  can  be  evaluated.  Using  the  boun¬ 
dary  values  for  at  the  end  points  it  is  found  that 


Gt(o)  =  o  =  A 

GtU)3  O  -  b  cos  \i  \  £  \J| 


(287) 


(288) 


Substituting  the  above  values  back  into  equation  (286) 
becomes 


G<x,f)  = 


B^nXx  X<f 

Dc es\x-  \>cos\[  sin k>3 

*0  x4 


(289) 


Because  the  Green's  function  is  continuous  across  the 
entire  interval  G,($)  «  G^s)  .  Using  this  condition  in 
equation  (289)  above  will  further  reduce  the  unknown 
constants 


B  S»|)  X§  s  t>  COS  X§  -  b  COS  yg  sm  >  P 

sminr  * 

B  =  b  cosX§  -  b  cos  Xl 

svoXl  sm  Xl 


(290) 


(291) 
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and  expression  (289  )  becomes 


WO.  VI*  \  <  r  l  J  ^7  /  WV  V  V1U  v  V 

0 [—41  -  — 5-^rl  *»>  **f 

(  *•*  X§  T^TvT  J 


t)  5 Cos  Xx  -  cos  Xj.  s»o 

(  ^nT 


} 


*>/ 


The  last  property  of  the  Green's  function  is  that  at 
the  discontinuity  Xr§  ,the  derivative  of  GO(|f)  has  the 
magnitude  fa)  or  -  |  by  expressing  equation  (282)  in 
self-adjoint  form.  The  discontinuity  becomes 


<>*  T*  "  1 

Using  the  terms  of  expression  (292 )  and  substituting 
into  equation  (293)  the  constant  P  is  found  to  be 


^  _  s\n  Xs 

X 

This  value  of  P  can  now  be  substituted  back  into  equa 
tion  (292)  and  after  several  algebraic  manipulations 
becomes 


G<x,o 


5,n  X(i-f) 

XsmXJ( 


XS 

X*m>4 


Xti-x) 


X4f 


(292) 


(293) 


(294) 


(295) 
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Appendix  C 


Determining  the  Green's  Function 
for  Case  III 

The  differential  equation  for  the  one-dimensional 
Cauchy-type  equation  is 

-2m  =  x4 

In  order  to  determine  the  Green's  function  for  a 
differential  operator  the  operator  must  first  be 
expressed  in  self-adjoint  form,  that  is,  it  must  ap¬ 
pear  as 


<*) 


Equation  (29 6)  is  not  expressed  in  self-adjoint  form 

2 

but  can  be  by  dividing  through  by  X.  The  Green’s 
function  must  satisfy  the  differential  equation 


-  Scx-s) 


(296) 


(297) 


(298) 


Because  G^f)  can  be  considered  as  two  separate  functions, 
G,(X)  and  G,(x)  ,  the  Green's  function  will  be  of  the 
form 

G6*f)c  fG,(x)  *<s  (299) 

For  and  will  each  satisfy  the  homo¬ 

geneous  equation  of  expression  (298),  that  is 
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(300) 


§w  -  fi G-‘"  = 0 


«U* 


-  2L  GXU)  =  o 


(301) 


The  solution  of  equations  (300)  and  (301)  can  be  ob¬ 
tained  by  substituting 

G  “  (302) 

from  which  ttU-M  can  be  found.  Therefore  the  Green's 
function  takes  the  form 


Ax1  + 

Bx* 

x  <  5 

Cx-'  + 

Dx* 

x>5 

(303) 


Utilizing  the  fact  that  G'y,*)  is  always  zero  at  the  end 
points  over  the  interval  (1»3)  two  of  the  constants  in 
expression  (303)  can  be  evaluated,  it  is  found  that 

6,(»)s  O  s  A+B  •*.  B"-A  (304) 

G2(3)  -  O  =  |  (305) 

Substituting  the  values  for  B  and  D  into  equation  (303) 
the  Green's  function  becomes 

Hi-**) 


GfxiS)  = 


(306) 


(  - - 

U  21  J 


*>/ 


Next,  using  the  continuity  property  of  the  Green's 
function,  namely 
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(307) 


&,<E)  *  G,l f) 


it  is  found  that 


c  («--£) 

••  A' 

Substituting  expression  (309)  back  into  equation  (306) 
for  G,(x)  will  give 

f  {c(*"-£)/(r-f))(  *‘)  x  <  ? 


(308) 


(309) 


G(x,f ) = 


(310) 


C(X--  t) 

2  n  / 


x>$ 


The  constant  C  can  be  evaluated  by  using  the  discon¬ 
tinuity  property  in  the  first  derivative  of  the  Green's 
function,  i.e., 


.  4£,<o  -  -1 

«)  x  <*x 


(3  U) 


Differentiating  both  terms  in  expression  (310)  and  sub¬ 
stituting  into  equation  (311)  gives 


7g  '  3  f 


(312) 


The  value  of  C  can  now  be  substituted  back  into  ex¬ 
pression  (310)  for  the  Green's  function,  which  after 
simplifying  becomes 
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Appendix  D 


Analytic  Solution  for  the  Two-Dimensional 
Steady-State  Heat  Conduction  Problem 

The  two-dimensional  Poisson  equation  representing 
the  internal  temperature  within  some  enclosed  region 
is  given  explicitly  by  the  partial  differential  equation 

V*T*  =  -  Rx.y)  (314) 

•if 

where  T  £x,y)is  the  temperature  at  the  point  (x,^)  and 
F(*.y)  represents  the  internal  heat  generation  from 
some  source.  The  solution  of  equation  (314)  will  be 
derived  over  a  square  plate  with  the  Dirichlet  boundary 
values 

T(x.o)»T(o.;i)=Tfx,^»  O,  T(«,3)  =  T.  (315) 

where  "1^  is  a  constant  temperature  along  the  edge  X  =  ^  . 

Applying  the  superposition  principle  the  problem 
can  be  simplified  by  assuming  that  the  complete  solu- 
tion  for  can  be  represented  by  two  independent 

harmonic  functions.  The  first  of  these  satisfies  equa¬ 
tion  (314)  but  with  homogeneous  Dirichlet  boundary  con¬ 
ditions.  The  second  harmonic  function  satisfies  the 
Laplacian  of  (314)  using  the  inhomogeneous  boundary 
condition  at  X-CL.  The  complete  solution  is  then 
given  by  the  superposition  of  both  harmonic  functions, 
see  Figure  25 »  that  is 
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■*v-” " 

•-I.  •  •• 


y 

>  o 

r‘ 

o 

T*  =  Tl+Ta 

o 

o  o 

=  O 

T'  V'f  ■  -  fix,) 

o  x-< 

X. 

O  X«< 

i 

Fig.  25  Method  of  Harmonic  Solutions  for  the  Poisson 
Equation 


T*=  T,  t  Ti 

or 


(316) 


^  (317) 

and  therefore  T*  satisfies  the  original  Poisson 
equation  hut  is  composed  of  the  two  harmonic  functions 

T,  and  ”Tj 

Poisson's  equation,  with  the  homogeneous  boundary 
conditions,  can  be  solved  by  assuming  that  can 

be  represented  as  an  eign-function  of  T,  ,  that  is 


-  =  XT,  018) 

Substituting  this  value  for  will  give  the  par¬ 

tial  differential  equation 

£l<  +  £!■  =  XT,  (319) 

By  using  the  method  of  separation  of  variables  for  7)  , 

Tf=  X'Y  »  the  two  eign-f unctions  are  obtained 
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i  •  B  *  $ 


(320) 


X(x)  =  f\s\t\n£%  YM=  Bsm  wffv 

a  J  T 


with 

t1  ) 

Therefore  TJ  it  is  found 


(321) 


T,  =  An*,  Sl*»  £p  S\nrft£Tj 

Because  of  linear  superposition  an  infinite  series 
of  equation  (322)  will  also  satisfy  the  original  diff¬ 
erential  equation,  therefore  Tt  is  of  the  form 


(322) 


7]  -  £  $1  Anm  sin  n/TX  sm  mrr  y  (323) 

n  €L  -p  j 

Substituting  expression  (323)  for  £  into  equation 
(318)  the  following  expression  results 


ftxiy)  -  A  £  £  A**,  sun  OJL*  sm  wff  y 


Solving  for  the  coefficient  Axm  it  is  found  that 

a  k 

"  aV~X  ^  J"  Sln  *12L*  Sinmiry  JxJkj 

©o  ^ 

For  Rfc«j)s  |  ,  the  coefficient  becomes 

A  -  1UU _ 

~  nrr»TTH(ttlw\'+ 

and  T,  ,  from  expression  (323 )»  becomes 


(324) 


(325) 


r 

\ 


(326) 


(327) 


T*  ?  “  I  £  S'n  "«P  S'r‘  ^  ^ 


The  second  harmonic  function,  "T^  ,  satisfies  the 
Laplacian  of  equation  (314)  and  has  the  inhomogeneous 
boundary  condition  of  at  x  =  q  •  The  Laplacian  of 


is 


(329) 


^T»  *  +  ^T*  =  O  (328) 

Again,  using  the  method  of  separation  of  variables  for 
^ Z  ,  the  two  eignfunctions  are  obtained 

XU)  *  A  *tnh  OIL*  <md 

Therefore  /z  becomes 

£  *  X-y  =•  Bn  Sio t)  inrx  sm  MTH  (330) 

Tr  T 

Using  the  linear*  superposition  principle  equation  (330) 
becomes 


Z.  *  SL  Bn  stnVi  mr  *  sin  nw  h 

n  T  “T 

The  inhomogeneous  boundary  condition  at  Xs*Lcan  now 
be  applied  to  evaluate  8»,  ,  that  is 

T*  T0  =  £  B^^mVmreL  smnxr^ 

b 

and  is  found  to  be 


(33D 


(332) 
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(333) 


+  To 

n  rr  $int)  aira. 

~~r 


Substituting  this  value  for  Bn  into  equation  (331 )»  Ta 
is  found  to  be 


5  HTo  s»nh  ^ 

n  W  *»nh  mra 

fc> 


n'n  ^ 


(334) 


The  complete  solution  for  ]  is  now  given  by  equations 
(327)  and  (334) 


T  *  W£  « 


-  £  y  r$±  s»n 

^79  n  nr*  _ _  b 


***(<?*?■  +  \>zr?) 


»  MTo  ^  nn  x 

IT  ^  _  *> 


nfro- 

W 


.  s\n  httm 

b 


(335) 
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Appendix  E 


Derivation  of  the  Cubic  Spline  Function 

The  spline  function,  ,  is  a  polynominal  of 

degree  three  or  less  on  the  subinterval 

I*.  *  t  »«■„«.] 

and  will  satisfy  the  conditions  that 

and  S‘(x;)*f(X;)  (t=k-»,lO 

Equation  (337)  imposes  the  conditions  that  the  spline 
function  have  the  same  values  as  fio  and  its  first 
derivative  at  each  node  location.  This  will  allow 
S(x)  to  be  substituted  for  f(x)  on  each  interval, 

To  insure  a  smooth  fit  of  adjacent  splines  across  the 
entire  interval  it  will  be  required  that 

=  SkJ*k.) 

and 

&V)  =  SJv) 

Expressions  (338)  and  (339)  assure  a  smoth  fit  of  ad¬ 
joining  splines  at  each  interior  node,  X#  .  Since 
SKU)  is  a  cubic  polynominal,  S^x)  is  a  linear  func¬ 
tion  of  X  on  the  interval  X  ^  .  The  second  deriva¬ 

tive  of  the  spline  will  appear  as 


(336) 


(337) 


(338) 


(339) 
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(340) 


s^ft)  -  ^(x*-*)  +  Wk(x-  XK.,) 

where 

Mk=  5^Xk)  (341) 

Expression  (340)  can  be  obtained,  from  Figure  26,  by- 

solving  for  'V’O  using  the  linear  slope,  ^*V  - , 

^  ic  *  ^ 

across  the  interval  (  y.  v  \ .  k_l 

V  XK-«  >  Xk) 


Integrate  equation  (340)  to  find* 

=  -MK.,(X|;-Xy-  t  MK(WK.y  +  C|  t»i 

2  2  Vj^ 

where  is  a  constant  from  the  integration  and  • 

Integrating  again,  the  spline  function  is  found* 

f  C4x  +-  Cfc  (343) 

4  hir  6  hK 
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Let  *jK  *  *  and  using  condition  (337).  <jK=  S^(vk). 

Therefore 

5  +  c,Xw  t  Cz  (3^) 

(p 

and 

Sk(xk)=  £UV  +  c,xK  +  Cl  (345) 

Solving  for  C%  and  C ^  s 

C.  =  (bK-3K.,)  k-,)(^k/t)  046) 

and  ^  ^ 

C*.  =  ~  x*-«bO  ~  (**mIm  -XicMuK^I^  (347) 

g  ~  -  - - — ^ 

Substituting  expressions  (346)  and  (3^7)  back  into 

equation  (3^3)  the  cubic  spline  becomes 
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Appendix  F 


Integration,  Using  the  Cubic  Spline  Function 

The  integral  representation  for  the  temperature 
along  the  line  is  given  by  equation  (264)  as 

=  I  GCx(§)  P(S,yn)4§  +- 

o  * 

The  one-dimensional  Green's  function,  G(x,0  is  given 


(349) 


by  the  expression 


G(x,$)  = 


x(a-s) 

a 

f  (q-  x) 


x  <  5 


*  >f 


(350) 


Substituting  in  ,  over  the  appropriate  limits, 

(349 )  becomes 

& 

+  l  P(5'9«)<is  +  (351) 

The  value  of  the  temperature  at  the  node  Xj^  is  found 
by  substituting  X  =  into  (351)  above  and  is  expressed 
as 

XK 


a 

J  +  TsJLk 


(352) 
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Since  a.  and  X^  are  constant,  and  because  ptf.9n)  is 
only  a  function  of  the  variable  §  for  constant  , 
equation  (352)  simplifies  to 

T(*k,i,n)  -  (~^)  J§  +  *,cjp(s)<l§ 

fp'nis  +  T^Xk  (353) 

*k  <*■ 

Now,  will  approximated  by  the  cubic  spline  rep¬ 

resentation  given  by  equation  (258),  substituted  back 
into  the  integral,  equation  (353),  and  integrated.  The 
results  after  much  algebraic  manipulation  yield 


Replacing  the  integrals,  given  in  equation  (353)  by 
their  equivalent  quantities  above,  equation  (353)  be¬ 
comes 


( 


(356) 
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